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ABSTRACT 


Unsteady aerodynamic modeling techniques are developed and applied 
to the study of active control of elastic vehicles. The problem of 
active control of a super-critical flutter mode poses a definite design 
goal--stability, and is treated in detail in this thesis. 

The transfer functions relating the arbitrary airfoil motions to 
the airloads are derived from the Laplace transforms of the linearized 
airload expressions for incompressible two-dimensional flow. The trans- 
fer tunction relating the motions to the circulatory part of these 
loads is recognized as the Theodorsen function extended to complex 
values of reduced frequency, and is termed the generalized Theodorsen 
function. A brief critique of previous attempts to generalize the 
Theodorsen function is given. Inversion of the Laplace transforms yields 
exact transient airloads and airfoil motions. Exact root loci of aero- 
elastic modes are calculated, providing quantitative information regard- 
ing subcritical and supercritical flutter conditions. 

The technique of generalizing simple harmonic airload calciilations 
to complex values of reduced frequency is extended to compressible flow 
regimes. It is conjectured that computer programs which calculate air- 
loads for oscillatory motions can be generalized in a fairly straight- 
forward manner to calculate airloads due to arbitrary notions. This 
is accomplished for the two-dimensional supersonic case. 

The ability to calculate airloads for complex values of reduced 
frequency allows approximate techniques of calculating these loads to 
be evaluated. Matrix Fade approximants of airloads for two-dimensional 
airfoils are evaluated in this manner. 

The exact airfoil motions contain portions associated with rational 
transforms and portions associated with nonrational transforms. The - 
oscillatory response characteristic of a fluttering airfoil is associated 


i 
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with lh(^ rational portion and a Uu‘ort‘in is provocl rr^ardinK ih^' cmmi- 
striic t ion oi a uniquo i inito-dimcnsional , I incar, consian t -c'oc / i ic i<*n t — 
modt'l t)l this portion ot the system. This rational modeJ does no( 
require' state aii};’men tation to model unstcaidy aerodynamic eiteerts and may 
be used to desiRM active aeroelaotic control systems. 

The rational .model and Fade model are used to design ilutter suppres- 
sion systems for airfoils in incompressible and supersonic flows using 
the optimal regulator design technique. Both techniques are shown to 
produce valid flutter mode control designs. 
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Chaptex’ I 
INTRODUCTION 


During the past decade, an aircraft design philosophy has emerged which 
attempts to gain performance improvements by means of an interactive de- 
sign process involving structural dynamics, aerodynamics, and control 
systems. This philosophy involves the use of active control systems to 
achieve aerodynamic and/or structural designs which have better perform- 
ance, stability, or economy than can be achieved with conventional passive 
techniques. Many of these concepts have been implemented in the B-52 
load alleviation and mode stabilization (LAMS) [Ref. 1], and control con- 
figured vehicle (CCV) [Ref, 2] programs. The concepts used in this control 
configured vehicle philosophy include: augmented rigid body stability,, 

maneuver load control, ride control, fatigue reduction, gust allevia- 
tion, and flutter mode control. The first five items involve the static 
and dynamic performance of the flexible aircraft. The design goal of these 
items is typified, by the C-5A active load distribution control system [3] 
which was designed to reduce the wing root bending moments experienced 
by the aircraft and thus increase its service life. 

The last item, flutter mode control, is fundamentally different from 
the others in that the structural stability of the flexible vehicle is 
involved. While loss of the former items would result in degraded per- 
formance or a shorter vehicle life, loss of a flutter mode control system 
at a supercritical flutter condition would usually result in loss of the 
aircraft. Although the risk is high, the potential performance gains are 
correspondingly high and flutter mode control systems can be designed to 
reduce the structural weight of a vehicle or to increase the flight enve- 
lope of the vehicle by expanding flutter speed placards. Roger and 
Hodges [4] describe the flutter mode control system implemented for the 
B-52 CCV program and successfully flight tested, while Sandford et al. 

[5], document a system installed on a wind tunnel model. 
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rt^quii’c'd art} common to alJ of tlu*si‘ CCV pro- 


'IMio analysis toclmiciuo^i 
nrams and Involve tlu' study of unsti'ady aerotl ynaiuicsi inj^ arbitrary motions, 
strutMural dynamics due* to unsteady loading, and aorori ynamlci loading 
caused by control surl'acn? motion. The design of flutt(}r mode control 
systems places severe demands U|)on the analyst; the |)rimary design goal, 
is structuiml stability. Hence, this dissertation focuse*s attention upon 
techniques of analyzing flutter mode control systems. Of course, tJie 
techniques will also be applicable to the other CCV concej^ts. 


A. SURVEY OF LITERATURE 

The finite element method of structural dynamics is well developed 
[6] and will be assumed as the basts of the analysis of aircraft struc- 
tures, The infinite dimensional spaces required to describe solutions 
are reduced to finite dimensional spaces by the familiar technique of 
truncated normal modes [7], 

The study of unsteady aerodynamics has progressed along two direc- 
tions : 

(1) The calculation of the indicial loading due to impulsive motion; 

(2) The calculation of the loads due to simple harmonic oscillations 
of the wing or section. 

The former area was first investigated by Wagner [8] for two-dimensional 
incompressible flow. R, T, Jones [9] and Lomax et.al, [10], continued 
this line of investigation. A method for calculating the loads due to 
simple harmonic oscillations of a wing section in incompressible flow 
was first given by Theodorsen [11]* The corresponding solution fqr sub- 
sonic flow was given by Timman and Van deVooreu [12], and for supersonic 
flow by Garrick and Rubinow (13j, Methods for calculating the loads on 
throc-dimens ional wings due to oscillations of assumed or norinal mode 
shapes have derived from Possio^s integral equation |14], Techniques of 
treating the singularit ies of the kernel function ami obtaining solutions 
were given by Watkins et al . [1^1, have been extended to wings with 

control surfaces by Howe, et al,, [16]. Another calculation process, 


analogous to the f Lnito-olomont methocJ of structures, is the douhlet- 
lattico technique of Albano and Roddon [17], 

The prevalence of aerodynamic analysis techniques based upon the 
assumption of simple harmonic motions is undoubtedly due to the success 
of the theory in predicting flutter boundaries, Theodorsen and Garrick j 18] 
and Smilg and Wasserman [19], are representative of the methods traditionall 
used in the calculation of flutter boundaries. The latter reference intro- 
duced the concept of artificial structural damping. 

Attempts to extend Theodorsen’ s theory to deal with arbitrary motions 
(e,g., converging or diverging oscillations) were made by W.P, Jones [20], 
and by Luke and Dengler [21], Jones concluded that Theodorsen* s solution 
could be extended to diverging (unstable) oscillations but not to converg- 
ing (stable) oscillations while Luke and Dengler *s attempt to extend 
Theodorsen’s solution to stable motions was rejected in a series of articles 
[22] - [26]. 

The inability of U-g flutter analysis and oscillatory aerodynamics 
to give quantitative information regarding stable, subcritical flutter 
conditions [Richardson, 27],[Hassig, 28], led to methods of approximating 
this behavior based. upon convolution techniques, R.T, Jones [29] indicated 
the method of exponential approximation of Wagner’s indicial loading func- 
tion and used the convolution integral to obtain results for- arbitrary 
motions, Jones’ work was followed by Goland and Luke [30], Baird and 
Kelley [31], and Dugundji [32], Recently, Vepa [33, 34] applied the tech- 
nique of Fade approximation of oscillatory loads to derive expressions for 
loads due to arbitrary motion. Also, Morino [35, 36] has developed a new 
formulation based on the Green function solution of the governing partial 
differential equation which is valid for arbitrary motions. 

Whereas the ability to calculate airloads for arbitrary motions is 
of interest to the aeroelast ician for the insight gained concerning the 
approach to flutter, it is a necessity to the controls engineer who desires 
to design a flutter mode control system. The application of the design 
techniques of modern control theory requires that the plant to be controlled 
be described a mathematical model, preferably by linear, constant- 


roef r itncMil , orillnary d i (’ f\»ron tial equal: ioa^-i , Approx Imat i.un l:oci}uii quc’^i 
ba.sed upon exj^omnit ia 1 approxlma t tons to tndicial rc*sponse func Lions or 
Pade approxifiiants leas* natuiuuiy to such models In which the unstcaady 
aerodynamic el’tecL.s are simulated by augmented state variables. The li-52 
cev flutter mode control system was desigiuxJ using this type of mr)de] \ 2 \ 
and utilized th(' frt>quoncy domain control synthesis method. 

Optimal control theory is a well developed methodology for the syn- 
thesis of control laws to minimize a suitable performance index (Bryson 
and Ho, 371, Designs of flutter mode control laws using augmented state 
methods to represent the unsteady aerodynamics and implementing the 
optimal regulator solution are described by Turner [38] and Dressier [39], 
A program designed to study the active control of flexible aircraft which 
incorporates Morino’s aerodynamic theory is described by Noll and 
Morino [40] . Hov^cver, it has yet to be applied to a flutter mode control 
problem, 

A different approach was taken by Nissim [41] who developed a flutter 
suppression scheme based upon the concept of aerodynamic energy. A wind 
tunnel program testing Nissim^ s design technique is described by Snndford, 
et a 1 . [5 ] , 

Flutter mode control system designs are actually problems in distrib- 
uted parameter system theory, Wang and Tung [42] surveyed the field and 
references [43] - [48] typify the results of the theory. Sung and YU [49] 
present a formulation within which the flutter control problem can be 
treated, while Wang [50] presents a technique of stabilizing a system with 
a finite number of unstable modes which resembles the flutter problem. 


B. THESIS OUTLINE 

Chapter II presents the equations of motion of the typical section 
treated in this thesis and derives the generalized Theodorsen function for 
arbitrary airfoil motions. The Laplace inversion integral is used to de- 
rive loads due to transient motions and generalized unsteady aerodynamic 
loads are studied in comi)ress ible flow. 
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In Chapter III, the gcnorall /.eel loarln dovelopod In Chapter II are 
Incorporated into the equations of motion, and the locus of roots of the 
aeroelastic system is determined. The Laplace inversion integral is used 
to calculate exact airfoil motions due to flap command inputs, 

Chapter IV treats the problem of approximation of unsteady aerodynamic 
loads. R.T. Jones’ approximation of the Theodorsen function, and Vepa’s 
matrix Pad<< appi'oximants of compressible loads, arc compared to the exact 
solutions for aribtrary motions. 

The active control of aeroelastic systems is treated in Chapter V. 
Controllability and observability of such systems are investigated and the 
aerodynamic energy design technique is studied. The "rational model" is 
presented and compared to the Pad<? model. The models are used to design 
optimal regulator solutions to the flutter mode control problem. 

Chapter VI presents the conclusions of this thesis and recommendations 
for future research. 


C. SUMMARY OF COOTRIBUTIONS 

1. The problem of generalized aerodynamic loads due to arbitrary air- 
foil motions is investigated. The generalized Theodorsen function for in- 
compressible flow is derived using Laplace transform techniques. The same 
technique is applied to compressible unsteady airload calculations and 
results are presented for the case of two-dimensional supersonic flow. 

Exact root loci of aeroelastic modes are calculated and examples of exact 
transient responses, due to stable motions are given. 

2. The transient motions contain portions associated with rational 
transforms and portions associated with nonrational transforms. It is 
shown that the oscillatory motions typifying flutter phenomena are due 
entirely to the rational portion of the response. 

3. The generalized aerodynamic loads aie used to evaluate approxi- 
mate techniques for calculating these loads. It is shown that exponential 
approximations of indicial loading functions and matrix Pade approximants 
of oscillatory airloads provide valid models of unsteady airloads for 


i 


values ol complex reduced frequency near tlie J.uiaKlnary axis. 

4, The goiiorall/.ed Theodorsen function Is usimI to study static 
diverKcnce of typical S(^ctions. It is shown that static divergence 
corresponds to the emergence of a real iu)sitive pole of the system trans- 
fer function and occurs, In addition, to the original structural poles. 

5, Flutter mode control systems are investigated, The contrc.lla- 
bility and observability of airfoils Is studied. A theorem is given 
concerning the ability to construct a unique, linear model of the 
rational portion of the acroclastic system which does not require state 
augmentation. The resulting rational model and the Fade model arc used 
to design flutter mode control systems. 
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Chapter II 


UNSTEADY AEUOUYNAMIC MODELING 


A, TYPICAL SECTION EQUATIONS OF MOTION 

The typical section which. will be analyzed is shown in Fig, II-l, It 

has leading- and trailing-edge control surfaces which are aerodynamically 

unbalanced (hinge lines at leading/trailing edges), simplifying the des- 

t 

cription of the aerodynamics. Linear and torsional springs and k^) 

at the section elastic axis restrain motion in the plunge (h) , and pitch 

(a) degrees of freedom, while torsional restraining springs (k„ and k ) 

P T 

restrain control surface. deflections . All linear coordinates (x, z, h) 
have been nondimensionalized by the semi-chord, b. The equations of motion 
are derived in App. A following the conventions of Theodorsen [11], and 
Theodorsen and Garrick [51, 18] as 


where the subscript 
tural origin, and 


M X 
s- 


s 


= -K X - B X + ■ - L + Gu (2.1) 

S‘*~ S“" . ^ 

m b 
s 

indicates that the matrix operators are of struc— 


M 

s 


a 


r 


a 




a 


[rp+Xp(c-a)] 


[rp+Xp(c-a) ] 




I Jc [x (d-a)-r ] 

L r r ' r 


The matrices M , K ^ B , and G arc 
s s 


On an aerodynamically balanced control surface, the hinge line is some 
distance away from the leading/trailing edge such that the aerodynamic 
pressure distribution may be used to advantage in reducing the hinge 
moments developed by surface motion. 


- 7 - 







riio us(> of Ltnwivni suivt'aco sjunuK and thimpiiH* constants to a])|)roxima t o 
i rr('Vtn\s lb lo position contx’ol systems is discussc»d in Ap)i, A, Kfiuatlnii 
(2,1) dcsscribes a foui’ dof.»ree-of-f recdom (DOF) model. Two and three DOF 
models may be obt ained -JCx’om (2,1) by deletlnf^- nppxnxpriate rows and columns 
of the- matrices and vectors. 

The specification of the aerodynamic load vector, L, completes the 
system description and is the subject of the remainder of this section. 


JB. UNSTEADY AERODYNAMICS 

The development of the linearized, small disturbance partial differ- 
ential equation for unsteady aerodynamic loads is presented in numerous 
textbooks and the present at ion --of Bisplinghof f , et al. [7], will be 
followed. The exact, nonlinear, unsteady flow partial differential equa- 
tion satisfied by the velocity potential is 



and the flow velocity is given by 


( 2 . 2 ) 


q = V0 (2,3) 

2 

where the V and V operators imply the use of dimensional coordinates 

X* = bx; y* = by, z*^ = bz . (2,4) 

The flow velocity is related to the pressure through Kelvin’s equation 
or the unsteady Bernoulli equation 

I? -h (q . V)q = - Vp . (2.5) 

Equations (2,1) through (2.3) are linearized by assuming that the fluid 
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volouity vectoi', q, varies only slightly from the f roo-sti’eain velocity, 
U. A flisturbancc velocity potential is defined such that 


<}) (j)' + Ux” 


where the disturbance velocity com ponent s 


u' = u - U = 


S(|) ’ ^ » 

; V = , w = 

Sx* dy* 


are assumed to be small compared to U. Then the linearized partial 
differential equation for unsteady, compressible flow is 


2*. . 1. ^±1 . m A'- _ 

a! a Sx*at ax*^ 


(2.6) 


subject to the boundary conditions 




w 


az* 

U ^ n u 


-b s X* s: b 


(2.7) 


az* az* 

* L L 

" " ^ ' 


■b X* s: b 


( 2 . 8 ) 


where z*(x* t) and z*(x* t) describe the location of the upper and lower 
surfaces of the section as shown in Fig. II-2, The linearized versic'i of 
(2.5) gives the pressure coefficient 


c 

P 


P-Pe, 

ip 


2 ac})' 2 a(i)' 

at U dx'^ 


(2.9) 


yielding the pressures on the" top and bottom surfaces of the airfoil as 
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FIG. II-2 CROSS SECTION OF A THIN AIRFOIL 






P 


P 


u 


L 


Pco 


p«u — 0^ t) 
c)x* 

P09U ^ 0; t) 


— ct)' (x*, .0- , t) (2.10) 
" St 

P„ 4)’(x^ O', t) . (2.11) 


Since the governing differential equation, (2*6) is linear, the solu- 
tion may be constructed as a super-position of elementary solutions* The 
airfoil profile may be separated into a portion representing thickness, 
z* and a portion representing angle-of-attack and camber, 

t n 


z 


* 

U 




(2.12a) 

(2.12b) 


The thickness distribution, z*, represents a symmetrical airfoil at 
zero incidence and, by syimnotry, can provide no lift or pitching moment. 
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The distribution, roprosonts a camburod, aovg thickness inclined 

mean line which produces the lift and pitching moments acting on the air- 
foil, This distribution may bo further separated into a steady portion 
containing the aix'foil camber and a nonsteady, mean— line portion nominally 
at zero angle-of-attach. It is the latter, 'flat-plate' airfoil which is 
the starting point for linearized, unsteady aerodynamic theory. Hence- 
forth, z*(x* t) refers to this flat-plate airfoil and <1) will be the 
velocity potential satisfying (2,6) subject to the boundary condition 



s w-><- 


dz* 

+ U — f 

St dx* 


(2,13) 


The flow prescribed by this boundary condition is antisymmetrical 
with respect to the x-y pia^.c, as described in Bisplinghoff [52], and 
the perturbation pressures at corresponding points on the top and bottom 
satisfy ^ ^ Thus the pressure difference act- 

ing on the airfoil, positive for downward loading, is 

p(x*,t) = Py - Pl = -2p„ u <t(x*,0^,t) - 2p^ — <t)(x^O;t). ( 2 . 14 ) 

ox* c)t 

The aerodynamic loads acting on the airfoil are determined by integrating 


this pressure difference 

over 

appropriate portions of the 

airfoil , 







P 

L 

p(x*, t)dx* 


(2.15) 

a 

b 

c 

(x* - ab)p(x*, t)dx^‘ 



M 

u 


(2.16) 



f 

(x* - cb)p(x*, t)dx* 


(2.17) 


■'cb 





c 

(db - x*)p(x*, t)dx* , 


(2.18) 


^-b 
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TIio method solution of (2.d) dopcnds upon tlio aerodynamic regime 
under investigation. In incom])ross t])le flow, M (J, and Liu* equation 
reduces to Laplace’s ociuation 


= 0 (2.19) 

which is an elliptic partial differential equation. In subsonic and super- 
sonic flows, the equation becomes one of hyperbolic type. The solution of 
the partial differential equation has traditionally been simplified by 
assuming that the airfoil is undergoing simple harmonic oscillations in 
the various degrees of freedom, thus removing one of the independent 
variables, t. Further simplification of (2.6) results if two-dimensional 
flow is assumed, making the equations .i.ndependent of the span- wise coord- 
inate, y. 


£. TWO-DIMENSIONAL, INCOMPRESSIBLE UNSTEADY 
AERODYNAMICS FOR SIMPLE HARMONIC MOTIONS 

A solution of (2.6) was first obtained for the case of two-dimensional 
airfoils undergoing simple harmonic oscillations in incompressible flow, 
Theodorsen [11] was the first to publish the complete solution, although 
many other authors obtained similar results independently during same 

period. B ispLinghof f e t al , [52], and Garrick [53,. 54] present summaries 
of the various authors and techniques. Appendix B contains a summary of 
Theodorsen’ s derivation as presented in Ref, 52, The solution consists 
of a superposition of flows due to a source— sink distribution, a bound 
vortex distribution along the chord, and a wake vortex sheet distribution 
convicted downstream from the trai ling-edge . The Kutta condition of 
smooth flow at the trailing-edge is enforced by Eq, (B.16), 


Q 






(B.16 


y 
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iuiuaLLon H,i7, giviuK Ihv circulatory lift, is roprosont a t i vo of tlio 
integral equations involved in the unsteady loads 

1+0 ^ 

Pc " J r^(^.o)d^ . (13.17) 

1 * 

To proceed with the vSolution, Theodorsen assumed 

(1) The airfoil motion, w*(x*,t) consisted of simple harmonic 
oscillations (Eq, B.21;, producing the wake vortex distribution 
given by (B.22) ; 

(2) The motion had been sustained for an indefinitely long period, 
allowing the upper limits on the integrals in (B.16) through 
(B.18) to be replaced by <x> ^ 

It was then possible, using an integral representation of the modified 
Bessel function K^(s) (equivalent to Eq. B.28) to evaluate the ratios 
of integrals occurring in Eqs, B.24 and B.25 as 

K^(ik) 

T7~^ rilTT r~T * (b.32) 

qk) -!• ^(k) K^(ik) + K^(ik) 

The restriction on the use of the integral representation of K^(s), Re(s) 

> 0, is not mentioned in the early references dealing with the subject. 
This restriction, in connection with the assumed airfoil motions (B.21 and 
D.22), implies the oscillatory divergent motion and wake vortex distribu- 
tion shown in Fig, II-3a, The analysis so far presented would thus appear 
to be inapplicable to the convergent oscillations shown in Fig, II-3b, 

The fact that the theory agreed with experimental observations of flutter 
boundaries (c.g., Theodorsen and Garrick, Ref, 18) explains the acco]^tance 
of the theory for Re(s) r: 0 (purely simple harmonic oscillations), 
although the integrals upon which the theory is based are then divergent. 

The simple harmonic loads acting on the airfoil are given by (B,35) 


C(ik) = 


H 




-in- 


✓ 


and (B . 36) as 


IRFOIL OSCILLATIONS 




f 




2 2 
pb U 


{-u 


M l-lk(B +c(ik)RS,J 

nc nc 2 


(K lC(:i.k)ltB, )lx 
^ nc ^ 1. ) -o 


( 2 . 20 ) 


1 k L* 

whor(' _x(x,t') • Kquation (2.1) spociali/.od to simple' harmonic 

oscillations and with u 0 is thou 


(M - 7]M ) h i 

s nc 




• U) k 


( 2 . 21 ) 


^0 


where T| = l/^^b. 

For a section with n degrees of freedom, (2.21) is an nth order 
matrix equation which has a nontrivial solution only if the determinant 
of the matrix of coefficients is zero. For a given airfoil section, (2.21) 
is a function of u, co, and k, and the determinant of the coefficients 
yields a complex equation. A method of solving this equation is to 
assume values of \x and k (allowing the aerodynamic loads to be calculated) 
and factor the resulting real and imaginary equations, giving two sets of 
n values of oj. In general, a given value of ca will not be a factor 
of both equations, and the process is repeated for other values of k 
until a value is determined for which the real and imaginary parts of the 
determinant have a common factor, , the flutter frequency. The small- 
est value of U corresponding to a solution is called the flutter speed 
given by -- This method of solution, termed Theodorsen*s 

method, is described in Refs. 11 and 55. 


An alternative method of solution which is more commonly used is the 

U -- g method described in Refs. 19 and 52. The artificial structural 

2 

damping, g, is introduced by replacing the real quantity, (l/co) , with 
the complex factor 





( 2 , 22 ) 


For a given choice of u and k, (2,21) now represents a complex eigen-* 
value problem for the unknown, Z. With the eigenvalues, Z, determined, 
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tho correspoiuliag froquc'ucy, vol(3city, ami structural clampini^ arc dotor- 
mined by 


2 

I 

Rc(Z) 

(2.23) 





J. 

k 

(2.24) 

g 

II 

(2.25) 


The critical flutter point is determined by the values of u and k yield- 
ing a value of g equal to the assumed structural damping (usually zero). 

The U - g method is commonly used to perform flutter calculations 
for compressible flow in which case the Mach number is an additional inde- 
pendent variable. In this case, the calculated flutter speed may not 

correspond to the density (altitude) and Mach number assumed in perform- 
ing the claculations and the analysis must be repeated at several Mach 
numbers so that a ^matched flutter point'' may be determined by crossplotting 
the results. 

It is obvious that a great deal of the calculation required to deter- 
mine a flutter point is of limited further value since the conditions 
corresponding to the intermediate solutions are unphysical. P'urther, the 
resulting flutter boundaries give quantitative results only for the case 
of neutral stability (simple harmonic oscillations). The information 
concerning the subcritical and supercritical flutter conditions is quali- 
tative at best. The cause of this situation is the assumption of simple 
harmonic motion in the calculation of the unsteady aerodynamic loads. 

Hence, an investigation of the possibility of -^calculating airloads for 
arbitrary motions is appropriate. 


D, THE GENERALIZED THEODQRSEN FUNCTION 

Attempts to generalize the Theodorsen function by evaluating C(ik) 
for complex values of k were made by W.P. Jones | 20] , and Luke and 
Dengler [21 J. Jones concluded that C(ik) could be generalized lor 
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(UvorR'lnK alrCoil motion:' (sou FIr. II-3n), but not for (lami)ocl converRliiR 
motions (FJr. II~:ib), Still, on tho IjuhIs of numoricnl calculations and 
clalminR that C(ik) could bo analytically cont;lnuod into tin' loft-haif 
plant:, Luke and DouRlor published tables of C(sb/U) for s ^ 0 i tw, 

0 < 0. Ilowt'vor, they did not offer a proof of this extension and in a 
series of re|)lles [Van do Vooron, 22 | , [Laltone, 23], IW,P. Jones, 241, 
[ChanRO, 2(5], their claim was re.lected. 

Earlier, Sears [56] liad used the technique of Laplace trap -formation to 
obtain new derivations of indicinl load-ing functions. Sears' presentation 
is essentially a derivation of the Rcneralizod Thcodorsen function although 
this is not discussed in Ref. 56. No mention is made [56] of criteria 
for the existence of integrals nor of the evaluation of C(ik) for complex 
values of k. 

The generalized Theodorsen function will be derived in a form closely 

following Sears [56]. Assume that the airfoil was undisturbed prior to 

t = 0 [w*(x*,t) =0, t < 0] and that the airfoil motion has endured for 

t ~ (x* - b)/U sec, producing a wake that extends from x* = b to x* ~ x* 

as shown in Fig, II-3. It is assumed that the airfoil motion w*(x*,t) 

a 

and the wake vortex distribution are Laplace transformable func- 
tions, Making the change of variables ~ (B.16) and (B.17) 

gives 



(2,26) 


(2.27) 


where c - (x* - b)/U. The change of variablCvS has the of fee L of making 
the wake vortex distribution function of a single variable 

Y"^^(t), LCjUations (2,26) and (2,27) are convolution integrals, and since 
the Laplace transform of the convolution of two functions equals the 
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p roduc l 


where 


oi’ tlu' t rap‘.) I’oriiis nf Ihp two Iuiu‘.tiou:> 


;; I'tiVl, Ki|S. (v.,2(l) aiul (2.2'n 


.(*, 1 q(i ) I 


— [v 


t I 2 lVU 




(2.2K) 


i'l.P^(t) ] - PU i'. 


t -I- b/u 
I- (2b/u)tJ 




n 00 

^ rr+ 2b/ u "* ^ ^ + 2b/u ' ^-st 


b s(b/u) 
U ® 




. S k(r> * \if[ 


rTT ^-(sb/«)t 


; Re(s) > 0 ; (2.30) 


t + b/u 


t + b/U 


+ (2b/u)tJ = ^ + (2b/u)t e dt 


00 t 


b sb/u 


U 


\y 


^1-77 dt 


i c"'”''" Kb(s) > 0 . ( 2 .:n) 


ORIGINAL PAGE IS 
OF POOR quality 



I 


111 ('vnimUitif; tUoHc oxiirosslons , tho ulurngo of varlnbJc.s 
aiiri (D.29) ami (15, 90) woi-c omployc'd . KliminatlaK f[y ] 
(2.29) ''' 


: (UL/b) + 1 

from (2.28) anil 


where 


and 


£l,P^(t)J = -2xcPbUC(S)r,[Q(t) J 



(2.32) 


(2.33) 


The Bessel functions in (2.33) are defined and analytic throughout 
the s-plane. except for a branch point at the origin and a branch cut along 
the negative real axis [Sect. 9.6, p. 374, Ref. 58], and by analytic 
continuation [57], C(s) is the unique operator relating Q(s) and L (s) 
throughout the s-plane (except along the branch cut). The principal 
branch of the Bessel function will be taken as -g < arg s s Jt and with 
the restriction on the real part of s removed, (2.33) defines the gen- 
eralized Theodorsen function. Setting s = ik recovers the Theodorsen 
function., (B.31). The remaining unsteady loads (M*^, and ) all 
involve the same ratio of integrals treated above, and the generalized 
Theodorsen function can be incorporated into the aerodynamic load ex- 
pressions by replacing C(ik) by C(s) in (B.33) and (B.36). 

For small values of |s|, K^(s) and K^(3) are readily calculated by 
their ascending power series expansions which are given in App. Jj. with 
b = re and C(s) = F + iC,. the real and imaginary parts of C(s) are 
plotted in Fig. II-4 which extends the figure given by Luke and Dengler 
[21] to 9 = + 60° and g = + 180°. The Theodorsen function is given by 
^he curves for 9 = 90° corresponding to the imaginary axis. As 
r -» 0, C(s) -4 1 and as r ^ c(s) ^ 0.5 independent of g. The maximum 
value of C remains in the range 0.2 < r < 0.25 Independent of o and 
increases monotonically as increases to 180"^, 
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Eiiuation (2,32) indicates l;Uat C(s) is to be ro^wr(\(u] as a froc|uoiu*y 
domain opei-ator or transit r function reiatini:^ Q(s) to L^ (s). Equ.ction 
(2,32) also proves that the Warner function and G(s)/s form a Laplace 
transform pair as implied by R,T, Jones [29] and Goland [59], and proven 
by Sears [56], For WagnerLs problem, Q(s) :: l/s anci 

. .2s£S£lii . 

^ s 

It is interesting to note that Scars’ development of this relation 
implicitly involved analytic continuation of C(s) through the deforma*- 
tion of the inversion integral contour into the left half-plane, although 
Sears does not comment upon this point. Although Sears states that his 
method is applicable to arbitrary airfoil motions, it seems that his 
intent was to perform such calculations via the convolution integral, 
using exponential approximations to the indicial load functions as shown 
by R,T. Jones [29], 

Equation. (2,32) indicates that the transforms of the aerodynamic 
loads will be multiple-value functions due to the branch point of C(s) 
at the origin. It is of interest to no^:e that Q(s) may contribute branch 
points also although this is not ^.he case for typical airfoil motions. 

Two additional derivations of the generalized Theodorsen function 
are available and are presented in App, D. The first derivation was 
given by W.P, Jones [20], while the second is based on the convolution 
integral. 

An outline of the calculation of transient unsteady aerodynamics, and 
a discussion of the difficulties in earlier interpretations of C(s), are 
offered in App, E, 

The Laplace transforms of the unsteady aerodynamic loads and the 
airfoil equations of motion for arbitrary motions are given by (2,20) 
and (2,21) v/ith ik replaced by s 

L(s) = +[B^^H-C(s)RS 2 ]s +[K^^+C(s)RS^]| X(s) , 
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(2.34) 


t 


IP* « 


X(s) = Gll(s) . (2.35) 


Equations (2,34) and (2.35) are matrix polynomial functions of s whose 
coefficient matrices contain the nonrat ional function C(s), (A rational 
functions of s is a function which can he expressed as a ratio of 
polynomials in s, [p. 60, Ref, 57],) In inverting these expressions, 
attention must be given to the branch cut of C(s) along the negative 
real axis. 


JE. INVERSION INTEGRAL FOR UNSTEADY AERODYNAMIC LOADS 

The time histories of the unsteady aerodynamic loads can be deter- 
mined from (2.34) using the Laplace inversion integral [57], To simplify 
the expressions, the unsteady lift, P(t), will be considered for two 
degree-of -freedom plunging and pitching motions. Equation (2,34) gives 
the transformed lift for this case as 


P(s) = 

where 

Q(s) = 


-jtPb 


2 U 2 1 

s h(s)+ (- s-as )a(s) -2«pbUC(s)Q(s) 


sbh(s) + Ua(s) + b(|-a)sa(s) . 


(2.36) 


(2.37) 


The inversion integral gives 


p(t) 


Oi+i“ 

f ^ 

I P(s)e^^ ds 


(2.38) 


where is to be chosen greater than all singularities of the integrand. 

The first expression in (2,36) is the noncirculatory lift, may 
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t 


be Inverted directly. The genernlized Thoodorsen function mny bo written 
in terms of the lift deficiency function 4>(s) introduced by Von Karman 
and Scars [60] as 


C(s) 

4»(s) 


1 - <))(s) 

K^(s) + K^(s) • 


(2.39) 


(2.40) 


Substituting in (2,38) 


P(t) = -jtPb^ ii + ^ a - aaj -2jtPbU hb + ua + b(i - a)o: 


- PbUi 


o^+i» 


Oj^-ico 


(2.41) 


<t>(i)Q(s)e®^ ds . 


The second term in (2,41) gives the "quasi-steady” lift P which results 

qs 

from ignoring the effect of the wake while the third term gives the effect 
of the wake. The integral may be -simplified by the deformation of the 
contour of integration [56] shown in Fig. II-5, The portions of the 
contour from a to b and from c to d lie above and below the branch 
cut of <J>(s) along the negative real axis thus making the integrand single- 
valued within the contour. The damped complex conjugate poles shown in 
the figure are representative of the singularities which may be intro^,. 
duced by Q(s). Sears treated the case of a step change in circulation 
[Q(s) = l/s], and proved that C(s)/s has no singularities within the 
contour given by N, N^,, N^, and the branch cut a-b, b-c, c-d. Thus 
the integrand .Is analytic at every point within the deformed contour and 
by Cauchy’s integral theorem [57], the integral around the contour is 
zero. If Q(s) is Laplace transformable, then the integrals along the 
semicircular arcs, and go to zero as the radius goes to infinity. 

The integrals along the cuts from--N and N to the poles cancel since 
the integrand is continuous along these paths, while the circular paths 

of infinitesimal radius around the poles give 2ni times the residues 
s t 

of [C(s)Q(s)e ] at the poles. Therefore, the integral- along the path 
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N is equal t:o 


)ds 


N 


j f(s)ds ^ { f(s)ds - { f(s)ds + { X(s) 

+ 2jti Res[f(s)J + 2jci Res[f(s)] 


(2,42) 


s=s. 


s=s* 


where f(s.) = 4 >(s)Q(s)e®^ , Along be, ds r: 0 ie^^d<) and the third 
integral on the right-hand side is 


(” ,(§5 

- -rr 


which approaches zero as q -* 0 if Q(s) ->00 no faster than l/s as 
s -♦ 0 which will be assumed henceforth. Along ab, s = re^’^ while along 
cd, s = re and 

J »(i)Q(s)e='a, = - j“[4.(f .l*)-,(££e-‘>')]<S(rc‘’')e-'‘dr 


'N 


' (2.43) 


+ 2ni|Res <■’( s)q( s)e®^j +Res |<)(s)Q(s)e®^J^_^^ 


Using the expressions 


K (re^^^) = K (r) + nil (r) 


K. 


^(re“^^) = -K^(r) + nil^(r) 


K (r)l.(r) + I (r)K (r) = l/r . 

o 1 01 


Sears [56] showed that the integral in (2,43) was 


K;'*/ * ’''(VH)T' 

o 


) «.rt , 

e dr 
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whtM’L' rt)/U is impliod as t:lu^ ar^uinc'iit. of Uu> Bcssc'I finuM, ions. Finally, 
the unsluady llTt is 


P( l) - P 1' P I- P 
' nc r nr 


(2.44) 


where 


P == -irPb 
nc 


H d . .o] 


(2.45) 


n , ■ 

-2rtPbU hb+Ua+b(i-a)d - ^ Res (<!>(i)Q( s)e®^) 

i=l ®“®i- 


(2.46) 


P = -2nPbU 
nr 


i 7 

•'o r 


, ijtv -rt 
Q(re )e 


(K - K -1- 
O i 




dr 


(2.47) 


and r = rb/U. and symbolize the rational and nonrational por- 

tions of the circulatory lift. The rational portion P^^ is comprised of 

the quasi-steady circulatory lift P , and a portion due to the residues 
s t 

of (t)Qe at the poles of Q(s), The summation in - P^^ is over all 

poles of Q(s). Typical airfoil motions result in rational functions for 

Q(s) which may be expanded by partial fractions into a sum of elementary 

transforms. The residues at the poles of the elementary transforms may 

then be calculated and used in (2.44), Table II-l lists several standard 

functions, Q(t), their transforms Q(s), and residue sums required 

in the evaluation of P (t). 

r 


Table II-l 


ELEMENTARY FUNCTIONS AND CORRESPONDING RESIDUES 


Q(t) 

0(.) 


6(t) 

1 

0 

1 

l/a 

0 

^-ot 

l/s-t-a 


Q wt 


sin wt - G cos cot] 

-at . 

« cos wt 


e ^^[(I-F)cos wt : sin wt ] 
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In evaluating the ro.siciue for c^(s) == l/(j?fa) the contour must bo in- 
dented at s " -a giving semicircular arcs of inf i ni teslmal radius. 

The integral exprOvSslon for (2.47)^ cannot be evaluated analyt- 

ically for typical airfoil motions but its integrand is a well behaved 
function and the Integral may be evaluated numerically. Figure II-6 
is a plot of the denominator of the integrand 




(2.48) 


As an example of the use of (2.44), P(t). is calculated for the case 

of a single DOF plunge motion (c^ pinned) with 


Q( t) = Ue ^*^sin u)t . 


The plunge motion yielding this function is 


h(t). = sin wt + w cos cot) 

b(cT +CJ ) L 


(2.49) 


and (2.44) through (2.47) and Table II-l give 




p(t') 


^ (iPU^)(2b) 


- 2;r|ie (5sin cot’-co cos cot*)-e (Fsin cot’4-Gcos cot’) 

- f [(K„-K, a-. 1 

•^O -* r[(-r-ta) +C0 J ) 


(2.50) 
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FIG. II-6 PLOT OF I VERSUS 




I 


whore? 0 - 0b/U and u> (ob/U, Tlio thrt>c oxprn.ssi ons on t;l»o right of 

(2.50) aro c > ^0 < r(.\spectlvoly, Figuro II-7 shows tiu? total 

■^nc •^nr 

and componont lift coefficients for 0 i. O. 2 IJ 2 corresponding to a 

damped oscillation with 0.7 damping ratio, and 0.2 rad natural fre- 

quency. Since Q(t) is continuous at t 0, the cii’culatory lift must 
start at zero which requires that the rational and the nonrational por- 
tions of the lift cancel at t = 0. Figure II-7 shows this to be the 

case, with c^(0) = c^ (0). The nonrational portion c , decays 

nc nr 

quickly from its starting value for small t' but decays slowly to zero 

for large t' and is a monotonic function of t'. 

A second example is a single DOF (degree of freedom) pitch-motion with 


0:(t') = 1 - cos wt’ . 


(2.51) 


For this motion 


Q(f) = 




Ujl-e cos cot* + ( 2 - 0)0 (dcos cot* + co sin cot’) 


and the resulting lift coefficient is 


c^ = 2 


n|-l+e I ) - t wt’ + - l)+ | pjsin wt* 

1 -5t* / f .-2 -2, _1 . ^ . \ 

2 e I ja(ui -a )-aJcos o?t* -(2ao)a + o)sin ct’ 1 




-1 


(| .i)(.?+5)4±’ 


/ - -^2 -2 
(-r+cr) H-u) 


-2n 

(-.j 

2 


- 

e dr 





i 


l'M.Kur(' 11.-8 shows tlic^ total and component lil’t coel’f lc::i cuts lor o m 

0. 2/^/2 and f\ ^ 0. l'’or this case tlie value ol Q(1) at t 0 is non/iTo, 

Indicating that tl\e circulatory lift does not start at /.ero. This is 

evident in Flu', 11-8 where it is seen that (0) f c^ «taj-ts 

at a slightly larger value than c^ .(0). Again, co decreases monoton- 

II r 

icnlly from its starting value and decays slowly to v-ero for largo t'. 

At t' -- dO, c^^_ has settled to within 1% of its final value, while 
has settled only to 16?o of its final value and contributes 5% of 
the total lif t . 

From these two examples it is clear that the nonrational portion of 
the loads will dominate at large t' for stable airfoil motions. Thus 
the asymptotic behavior of the loads is of interest. Sears [56] studied 
t)ie asymptotic behavior of the lift in Wagner's problem (step change in 
circulation) by using series expansions of the inversion integral integrand 
for small s since the behavior of P(t) for large t is determined by 
the behavior of P(s) for small s. The nonrational portion of the lift 
is given by the last term of (2.41) 


. O' -{-ICO 

^nr^^^ 1 . '*’(s)Q(s)e®^ds . 


(2.52) 




Using the ascending power series for and given in App. F 


t(s) = 


J- + •••](log I + 7^) + + ... 


1 + 


( - -2 - -3 

i — 1 V \ l ^ . Ay ^ ^ S 




(log 5 + - 1 ^) * |-(108 f + f-(iog I r_-i)t ... 


•Sdosf + r.) - S^iog| + I ^ 


- f + V I ' P'"-' I + 

+ log .'i] ^ ) . 


(2.53) 
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i 


r 


If Q(s) is of thp fona 


Q(«) 


, -m 

]) s j . . i I b s H b 
ni 1 o 


n s 
n 


a „ H 
n-1 


i- • • • {- 


a^ s 


III < n 


(2.54) 


thc'n for sninll ,s 


o(s) = + ••• I b^s I b^) i-(n^s i- ... + 


-ris2 


+ (ns + ... + n 5 ) 

1 n 


and the leading terms of <)'(s)Q(s) are 


♦ (i)Q(s) = -b^8(loe I + ”)'^) - (b^-b^a^)s^(log | + 


- b^i^( log I + + 


-n 


If Q(s) - s then the leading terms of <t>(s)Q(s) are 


(2.55) 


<t>( s)o( s) 


-S"*\iob|. V ) 


_n+2, s 
s ( log - 


ry 

e 


(2.56) 


Thus the asymptotic behavior of the loads is determined by operators of 
the form s [log(s/2) + Ypl'". Sears evaluated (2.52) using this oper- 
ator Integrated along a deformed- contour . Soars' development involves 
the questionable stop of utilizing expressions derived from the gamma 
function, r(n), evaluated at negative integer values (where P(n) is 
singular). Hence an alternative evaluation will be given which leads 
to the same results. Following Sears, the contour in (2.52) may be de- 
formed as shown in Fig. II-5. The asymptotic lift for stable motions is 
given by the nonrational portion (the integral in (2.43) with '!(s)Q(s) 
given by (2,55) or (2.56). The expression which results from the operator 
s^' log(s/2) + y l'" is 

ORIGIN AT. PAGE IS 
OF PouJt QU^^ nr 


( 


i 


i 






{ 


since the 
vanishes . 


while for 


since 



integral around the infinitesimal circle about the origin 
For m = 1 this reduces to 


“ 2^n I - in + dr 




(log •§ + in + T dr 

M e 


t' 3 ^ 1 


m = 2 the expression becomes 


2 ( 


2 ( 


- 1 )" ( 

n 


r -ft' 

log - e df 


00 

zii! f 


t' •'o 


e ^^(log - log 2t’ )df^ 


[r’(n) - log(2t* )r(n)] 


r'(n) 


= f 


.CO . -Tv 

-n-1 .. - . 1 , 

r^ (log r^)e dr^ . 

o 
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Thus, !ill;lu)UKh s" ^ | lon(;;/2) i yj"’ tioes not possc'ss :i I.nplncc transform 
for n ;> 1 , the; iiiv(M-so I.apiaco transform nivini by (2.52) can 1)0 ova I u- 
at('cl asymptotically as t' -» c», and s 0 In terms of this expression. 
'Hu* corren; pondoiK!^ is 


[loB I , t] - 




n ^ 1 ; 


(2.57) 


s"“^[log| + rj - lr’(n)-log(2t' )r(n)], n ^ 1 . 


(2.58) 


Sears evaluated the asymptotic lift using the expression given in (2.56) 
with n - 1 (Wagnei''s problem) and obtained 


c ~ 2 J-J- + (2f) 2 ^ \ 

nr ' t' t' ' 

Thus for step changes in Q(t'), the lift approaches its final value 
asymptotically as l/t ' . For the airfoil motion given by (2.49), the 

asymptotic lift is found from (2.54) with b = il)/(a^ + cj^) , b 0, and 

- / " 2 -2 ^ ^ 1 * 

” 2a/ (a + (a) ) giving 


nr 


-2jr 1 

1 3 
1 3 

1 

2 U) 

-2! 

7 +01 1 

Lt2 -2 _2 

[t* rf +a, 


■^3 


(r>(3) - 2 log(2t«) H- 


Thus stable airfoil motions for which Q(s) is of the form of (2.54) 

2 

approach final values asymptotically as l/t' . In Fig. II-7 c 

- ,-2 -2 " ’ 


(-oj/'iy + (jj )l/t' , while in Fig. II-8, c. ~ l/t’. 

*nr 


-^nr 


One further case of interest is that for which ()(s) - s“ with 
n s - 2. An airfoil motion of this type is jy = t'(n = -2). Tiien the 
leading term in (2.56) has an inverse Laplace transform given by 


■'[^= 


( log ^ + r ) 


= log t* 
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2n 1 ofT t ^ 


and tho asymptotic uonrational lifo coefficient in c. 

' f-nr 

while the rational lift coefficient is p^iven exactly by c^^^ - -2jtt* , 


K, GENERALIZED CQMPliESSIBLK AERODYNAMIC LOADS 


E, 1 General Eormulation 

When the flow is assumed to be compressible, tho Mach number M be- 
comes an. additional independent variable and the governing partial differ- 
ential equation, (2.6), is a hyperbolic equation. . Solutions [15], [12] 
have been obtained by assuming simple harmonic motion and making the 
substitution 4^(x, t), = which is equivalent to applying 

the Fourier integral transform [61] to the time variable of (2,6), In 
attempting to derive solutions for generalized motions, it is natural to 
apply the Laplace integral transformation. Defining 




oc 

c,z,s) = ^ 


<J)(x z, t)e dt . 


(2.6) and (2,13) become 


C> + (1) - — (1) - ^ -MO 

XX zz 2 ax XX 

a « 

GO 

= “ ^ <l)(x,z,0)- (t>^(x,z,0)- ~ q>^(x,z,0) ; 

a a « 


(2.59) 


X[w^(x,z,t) ] = <l>^j = sX[z^]-z^(x,z,0)+U ~X[z^], -bs:xs:b. 

’ z=o 


(2.60) 


(The variables x and z ire assumed to be dimensional throughout this 
section.) Let ft = 0’ + 4>" with 4>" a known function to be chosen such 
that 


2 


2Ms 




4,". - M d'' 

^xx zz X XX 

a 00 

00 

= - ft(x,z.0)- <t> (x,z,0) - — ft (x,z,0) , 

GO 


( 2 . 01 ) 


00 
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(I-*' 

z 


z-o 


(2.62) 


=5 -z (x 
a' 


z,0) 


The equation for is then 

2Ms „2^, 

~ 1’’ - - M (t» * 0 (2.63) 

a 

Gd CO 

4>;i = sX[z J + U 1^ X[zJ (2.64) 

Z=sO 

which is formally identical to the simple harmonic motion problem with the 
replacement of iw by s. 

Equation (2,63) is a homogeneous equation for whose boundary 

condition, (2.64), is linear with respect to Z[z ], while (2.61) is an 

a 

inhomogeneous equation linear with respect to the initial condition 

<t)(x, z, 0) and whose boundary condition (2.62) is linear with respect 

to the initial condition z (x, z, 0). Hence the transformed loads L(s) 

a 

due to airfoil motions X(s) may be written as a matrix equation 

L(s) = K’X(s) + K"x(0) . (2.65) 

It is interesting to note that for airfoil motions for which 

" is identically zero and the entire 

solution is given by (J)’. Also, since stability of a linear, system cannot 
be a function of initial conditions, the flutter problem is solely depend- 
ent upon (J) ' . 

The formal identity of the equations satisfied by 1> for simple 
harmonic motion and by for generalized motion implies that existing 

solutions of the simple harmonic motion problem may be applied directly 
to the determination of <>* by the replacement of iw by s. Thus 
the Mathieu function solution of Timman [12] can be generalized to provide 
solutions to (2.63) and (2,64). 


■!>' + <!'' 
XX zz 


- 39 - 


I 


I 


It is anticipated that the decomposition indicated by (2.65) occurs 
in solutions based on the acceleration potential, ilf, since tl? and 
satisfy the same partial differential equation. Also, the generalization 
of the above Laplace transform method to finite v/ings in three-dimensional 
flow offers no difficulties. Thus programs based upon kernel function 
techniques [15], [17], [62], [16], may also be modified to calculate the 
Laplace transforms of generalized aerodynamic loads. It must be empha- 
sized that the resulting transform is not the total solution, but corres- 
ponds to that portion of the solution which is linear in the transformed 
airfoil displacement modes. 

F, 2 Generalized Unsteady Supersonic Loads 

In the case of two-dimensional supersonic flow, Garrick. and Rubinow 
[13] obtained the solution for the simple harmonic loads using elementary 
solutions of (2.6) known as source pulses and gave the loads for the 
three degrees-of-f reedom: plunge, pitch, and trailing-edge control sur- 

face. Hassig [63] extended Garrick’s treatment to cover leading-edge 
control surfaces. The loads, due to arbitrary motion, which are linear 
in X(s) may be obtained from the expressions given by Garrick by the 
formal replacement of k by -is as shown in tie preceding section, 

(The resulting loads do not include those portions dependent upon the 
initial conditions of the motion.) The velocity potential of [13] on the 
upper surface of the airfoil becomes 





a o 



(x-?) 




( 2 , 66 ) 


with the airfoil lying between x == 0 and x =: 1. 

Alternatively, (2,64) may be derived directly from (2,63) and (2.64). 
following the procedure of Stewartson [64] summarized in Ref. 7 [pp, 364-367], 
Stewartsonts- procedure of Laplace transformation on x applied to (2.63) 
leads directly to (2.66) with the recognition of*^ 

^ z, s) 
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Z ws - a, 
X f X 1 


rj - ^o<v> 


where = sM/a(M -1), Garrick^s solution in terms of is recovered 

by use of the relations 


) ; -Tt < arg s ^ | 


• ~ -n'*! 

^ irt < arg s ^ it 


and noting that the Bessel function is single valued. Thus the 


above inverse transform is 


J l-i 


(M - 1) 


verifying (2.66) as the generalized velocity potential for supersonic 


two-dimens ional flow. 


Following the notation of Garrick [13], the axis of rotation is 
located at x = and the control surface leading edge is located at 
X = x^, Generalizing the expressions for the loads given by Garrick 
produces 


Pb(s) 


2 2 
8Pb U 


M^(s) 


- 1 S^[m (s)s^ + B (s)s + K (§) 


Mp(.) 


where 





M (k) 
a 


(q2-2x„r^) 




(f3~^^oP2) 




2 

— s 
3 3 


B (s) = 

a 


1 

2\ 


(q 


r^oq’ (2S■^’‘o<‘>l*^'■l'*^’‘^l) (^=2*“2'V’‘0>) 


2(P2-Vi> 


2s, 


K (s) 
a 


1 

— r 
2 1 


0 (q 




rVi^ 


The functional dependence of these matrices on i is meant to indicate 
that the parameters q^, r^, s^, and t^ given in Table II-2 are functions 

of i. All of these parameters may be derived from the ’Schwartz function 


f^(M, w) 


.1 — - 

I -ioou ^ , 

Jj- u)du 

- n 


(2.68) 


[where w -i (2 sM^/M^- 1 ) 1 by the recurrence relation 
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TJic Schwartz function (Ec|. 2.S8) i»s not expressible in terms of 
ciemcMVtary functions but it may be computed from the series 


-f (M,w) 
o 


-i") 

e 




(a)) I' iJ ,("0l 
ivid. ^ 

2^nl(2n I 1 ) 


given by von Box’bely [65] where the Bessel functions of complex argument 
are evaluated by their ascending power series given in App. F, Although 
transient time responses of the loads for specific airfoil motions, 

X(s), could be computed from (2.67) via the Laplace inversion integral, 
this has not been attempted. To perform this calculation would require 
knowledge of the singularities of the transformed loads which are not 
readily available. Note that the loads given by (2.67) do not involve 
a single nonrat ional transform such as C(s) for the M ~ 0 case. 
Fortunately, the exact transient time responses for indicial motions have 
been calculated by Chang [66] and Lomax, et al, [10]. The time responses 
of the loads for indicial motions at supersonic velocities are typified 
by discontinuous first deri’^atives and different functional dependence 
for various time zones. These facts indicate that calculation of trans- 
lent loads using inverse Laplace transformation would be laborious. Chang 
used the indicial response functions to calculate the simple harmonic 
loads from the convolution integral and noted [26] that arbitrary trans- 
ient motions could be treated in the same manner. The primary use of 
the transformed loads, (2.67), is for the investigation of airfoil 
stability — the flutter problem. 
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Chap tor III 


SOLUnON OF TIO^ AEROEMSTIC EQUATIONS OF MOTION 


Thu expressions for the loads, (2,34) or (2,67), may be substituted 
into the equations of motion, (2,1), giving 

G(s)X(s) ^ B(s) (3.1) 

with 

G(s) = [ms^ + Bs + K - Q(s) 

L s s s J 

and 

B(s) Gll(s) , 

Q(s) is derived from either (2.34) or (2.67). The primary goals of the 
analysis of (3.1) are the determination of the stability of the system 
and the calculation of transient motions. 

Since the airfoil is a linear system, its stability is determined 
by the homogeneous version of (3.1) 


G(s)X(s) = 0 , (3.2) 

Table III-l 

METHODS OF SOLUTION OF AEROELASTIC EQUATIONS 


Method 

Aeroelast ic Equation 

Solution 

Stability 

Criterion 

P 

{M + K - Q(s)}X(s) = 0 

s s 

s = "yi^+ik 

T 

k 

- -5 Q(ik)}x(il<) = 0 

(l+ig) 

g 

p-k 

{M + K - Q(lk)}x(E) = 0 

s s 

s “ yi^i ik 

7 

augmented 

states 

{M»s^ + K'-Q'(i)}x'(s) = 0 
s s 

s = ai Ig' 

f;=:rns ^ „ 


- 45 - 


OUIGINAI.?^ 





Noutrivial solutions arc glvtni by the zcroos ol tho cbaractor Is tit; <;quation 


i)ctLa(,s):i 0 (3. a) 

Which are the poles of the system. Table III-l Is drawn from Hassig [28] 
and summarlises the prevalent techaiquos of determining the /.eroes of this 
equation. The structural damping matrix has been eliminated for 

convenience. The matrix 0(s) involves structural, geometrical, and 
aerodynamic terms which influence the solution. The structux’al and 
geometrical terms are valid for arbitrary motions, while the motions for 
which the aerodynamic terms are valid depend upon the underlying theory. 

For instance , quasi-»steady aerodynamics (Q(s) = Q s + Q ) may be used 

Jl 

to analyze arbitrary motions for low frequency effects, while calculations 
of flutter boundaries commonly utilize aerodynamics which are valid only 
for simple harmonic motions. The p-method is intended for use in the 
former case, in which Q(s) is valid for arbitrary motions. If Q(s) 
is a rational function of s, (3,2) becomes a linear eigenvalue problem 
and solution by linear matrix techniques is possible,- Otherwise, the 
roots of the equation must be determined by iteration. The advantage of 
the p-method is that the exact roots and the degret' of stability of the 
system are determined, to the extent of the accuracy of Q(s), The 
stability criterion is that the real parts of the roots of the equation, 
yk, must be negative. 

The- k*“ met hod is the traditional U-g method which is used to determine 
the flutter boundary utilizing simple harmonic loads. Complex roots are 
obtained by introducing the artificial structural damping factor g, and 
a root of the equation represents a point on the flutter boundary if the 
corresponding value of g equals the assumed value of g. Disadvantages 
of the k-method are: (1) many solutions are required to obtain *’inatched- 

point” flutter boundaries; (2) for a given airspeed, several solutions with 
different frequencies may occur, leading to problems of sorting the roots, 
and (3) information obtained regarding subcritical and supercritical 
flutter conditions is only qualitative. Regarding the last point, 

Richardson [27] and Goland and Luke [30] give calculations illustrating 


( 


the ditforcnct\s between raten o£ change of damping at the flutter npeed 
calculated by the k-mothod and by more accurate mcthod.s. 

The p*-k method [67] attempts to Improve upon the k-method for su)>- 
critical and supercritical flutter conditions (i.o,, non-simple harmumic 
motions) by allowing the reduced frequency to be complex instead of intro- 
ducing the structural damping factor g. It assumes that if Q(ik) is 
calculated for oscillatory loads at s ^ ik, then the same loads will be 
good approximations to the true loads for s = yk + ik if y « 1 . 
Results given, by Hassig [28] confirm the usefulness of the method, 

The augmented state method is fundamentally different from the first 
three methods of Table III-l in that it attempts to model the unsteady 
aerodynamics with a rational transform. The primes on the matrices of 
the last row of Table III-l imply that the matrices have been modified 
to include the augmented states. The advantage thus achieved is that 
the resulting system may be analyzed by linear eigenvalue techniques. 

Note that the p-method has been used in the past with quasi-steady aero- 
dynamics (thus ignoring the effect of the wake) to maintain the rational 
form of the equations, while the k- and p-k methods have sacrificed 
the rational form to include the wake effects in more accurate oscilla- 
tory aerodynamics. The augmented state method is based upon R.T, Jones 
technique [29], [6], of exponential approximations to ind-icial loading 
functions, and wake effects may be approximated at the eKpense of the 
extra states. Since this form is wel uited to the needs of active 
control, it has found application in aircraft stability augmentation 
studies .-and is studied in detail in the nv9xt chapter. 

The differences in the stability criteria of the various methods 
may be delineated as follows: (1) the p and p-k methods determine a root 

at s = + ik; (2) the k-method determines a root at + i\^ = 

. 2 

(1/co )(l+ig); (3) the augmented states method finds a root at s =.• loo. 

The stability criterion for the p and p-k methods is given by the signs 
of Y = ^ where <J> is the angle between the ik axis and the 

root, wnilc the stability c riterion for the augmented state method is 
given by the sign of Q ~ •+• co^' corresponding to a characteristic 
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equation factor (s 2^ujs 4 * w ) [(s*ia) + w ] , For roots close to the 

ill) axis, C ^ T* stability criterion fox’ the k-methoU is deduced 

2 

from the characteristic equation (for a one DOF system) ms 4 . k (1+ig) = 
2 2 2 s s 

0, or s + igu) -f CO r: 0. Thus the, correspondcMice between g and Q 

2 

; 2^u)S igu) or g 2Q, The nonphysical nature of solutions based 

on the k-method for non/.ero values of g can bo seen by factoring- the 
characteristic equation: 


s 


-VFd - 

s 





^i(jj/2+<p/2) -T^^i(-]/2+(t>/2) 

* m ^ m ^ • 


The roots of the characteristic equation are not complex conjugate factors, 
emphasizing the unphysical nature of the solution. 

The remainder of this chapter extends the p-method of solution by 
using the generalized aerodynamic loading functions derived in the last 
chapter to study arbitrary airfoil motions. In the next section, the 
stability problem is 'studied and the last section studies arbitrary 
transient motions of airfoils. 


A, ROOT LOCI OF AEROELASTIC MODES 

In Ch. II, analytic expressions for Q(s) were given for two-dimen- 
sional incompressible and supersonic flow. The loads, Eqs, (2,34) and 
(2,67) are valid for arbitrary motions and give the exact airloads. 

They may be combined with the equations of motion as indicated in (3,1) 
and the stability of the aeroelastic systems studied via the p-method 
of Table III-l, Similar calculations are mentioned by Dengler, Goland, 
and Luke [25] in attempting to define their generalized Theodorsen func- 
tion but have evidently never been published. 


I 


Since the ioucl^; ai’O not rational funeticni.s o! , a eoni])uter program 
was flevolopod to numeric; a 1 1 y dotC‘rmin(‘ the roots ol Uie chni’ac tc'rist i <• 
equation, For the systems treated in this thesis, it was feaslljle to 
numf'rienliy expand the determinant in (3.1). A gradient search algoritlun 
was employed to locate the j/.eroes oi the detcnaninant which are Lhf> poles 
of the aeroelas-t*Lc system. The gradient was numerically determined by 
finite differences in the s-plane and the performance of the algorithm 
was quite satisfactory for the systems treated which included systems 
of four DOF (degrees of freedom), (eighth order). 

Table II 1-2 shows the behavior of the algorithm in converging to a 

root. Convergence is shown for a three DOF section and a four DOF section. 

The nondimensional velocity was U/ba) = 3.0, near the flutter speed of 

O' 

both sections and the search was started at s = 60 rad/sec. For the 
three DOF section, the algorithm locates the root to four significant 
figures in three iterations. Five iterations are required to achieve 
the same accuracy for the four DOF section. 


A~1 Incompressible Two-Dimensional Flow 

Table III-3 lists the structural and geometrical parameters of the 
three DOF system used in the following calculations. The frequency 
ratio while the natural frequency of the flap mode is 

three times the torsion mode frequency. 

The equations of motion, including the loads, were given by (2,35) 
and are repeated here: 


i 


(m -qM )s 
^ ^ s nc 





Gll(s) . 


(3.4) 


The loads contain the generalized Theordorsen function C(s) which is a 
function of s and U/b, Thus, with the nondimensional velocity speci- 
fied, the roots of the equations of motion may bo determined by iteration 
in the s-plane. Figure III-l shows the exact locus of roots of the three 

DOF system of Table III-3 as a function of U/bu) . The im^rtia coupling 
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oi' tho thr<n> modos causes tlio /.ovo airspeed uuturul frequeucles to be 
Hbiftod irom Lheir uncoupled values. As the airapeod increases, the bend- 
ing and torsion modes ai^e seen to approach each oLher in tho stable 
leit-halt plane with tlu‘ bonding branch becoming the flutter mode at 
U/b{j^ d.O, Tne flap mode remains stable throughout this speed range 
oven though - 0, 

Because the terms involved in the matrix of coefficients of (H,4) 
do not become infinite in the finite s-plane, the determinant has no poles 
and examination of the number of 360- phase changes of the determinant 
around a closed contour will directly indicate the number of zeroes 
within the contour [57, 61], This was accomplished for the section 

of Table III-3 at =3.0 by evaluating the determinant along the 

branch cut on the negative real axis, and along circles of radius 
0.0001 and 1000 rad/sec. Six 360° phase changes were obtained, accoiant- 
ing for the six known structural poles and it is concluded that these 
are the only isolated singularities of (3,4) at this airspeed. 

The root locus format is used for the presentation of results rather 
than the conventional U-g, U-co plots since the ability to calculate 
generalized aerodynamics makes this a more natural format. It avoids the 
numerical problems of root-sorting since the loci do not cross each other 
and it is required for active control design applications, 

A- 2 Supersonic Two-Dimensional Flow 

Table II 1-4 lists^ the parameters of the three DOF section used to 

illustrate the aeroelastic root loci in supersonic flow. The loads, 

(2,67), are functions of the generalized supersonic reduced frequency 
— — 2 2 

parameter, co z= -i2sM /(M -1), and the algorithm described above may 
again be used for the determination of the system poles. The equations 
of motion, including the supersonic loads, are 

+ j X(s) = Gll(s). 
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Table III-4 

THREE DEGREES OF FREEDOM SECTION PARA^^ETERS FOR SUPERSONIC FI.OW 



where T]' = - 1 . 

The locus of roots of this section are shown in Fig. III-2 as a 
function of Mach number. At M = 1.25, both the lowest frequency coupled- 
bending-torsion mode and the flap mode are unstable. It Is suspected that 
the flap mode is primarily a single DOF flutter mode [68]. As the Mach 
number increases, both of these modes become stable at M = 1.4. Above 
M 1,8, the remaining coupled bending— torsion mode flutters. Hence, 
for the mass ratio (a = 40, the range of stability for this section is 
1.4 s M ^ 1.8. 

The aerodynamic matrices M , B , and K are derived from terms 
composed of finite Integrals of exponentially weighted Bessel functions 
of integer order as shown by (2.68). Since these Bessel functions are 
single-valued analytic functions of s, there will be no branch points 
of (3.5) as in the incompressible case. However, a cursory review of 
supersonic indicial aerodynamics [e.g,, Lomax et a].,. Ref, 10], leads 
to the conclusion that (3,5) must nave more singularities than the six 
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structural polos, because those six polos cannot yield t\w complex indicia 1 
i'unctlons. SLneo the polos of the system arc the zeroes of the determi- 
nant of the matrix of coefficients of a search was made for additional 

zeroes of this func^tion, A circular contour of radius 1000 rad/sec centered 
at the orif^in yielded six 360° phase changes of the determinant, account- 
ing for only the six known structural poles. Further searching located 
the first additional zero as a complex conjugate pair at s = -1315 ^ 
i 1501, ovei ten times the frequency of the flutter mode. No other zeroes 
were located since the power series expansions used to evaluate the com- 
plex Bessel functions were numerically unstable at larger values of |s|. 
However, it is anticipated that an infinite sequence of additional zeroes 
of increasing modulus does exist, due to the oscillatory nature of the 
exponential weighting factor in the integrand of (2.68), and accurate 
transient response calculations would require the evaluation of a number 
of these zeroes of lowest modulus and their corresponding residues. For- 
tunately, the flutter problem can be studied by determining only the zeroes 
due to the str-uctural poles, as indicated in Fig. III-2, 

B. INVERSION INTEGRAL FOR ARBITRARY AIRFOIL MOTIONS 

Returning to the case of incompressible flow, it is possible to 
calculate exact transient motions from (3.4) using the Laplace inversion 
integral. With the substitutions 



»(s) GH(s) . 




I 

I 


K((ua*:ioa (3,4) becomos 

^(s)x(s) - (3,6) 

where Cf(s) and IB (s) are nXn and nxm matrices whose elements arc func- 
tions of s, and X(s) is an nxl state vector. The number of decrees 
of freedom of the airfoil is n while the number of control inputs is m„ 
If Det[a(s)] /: 0, [flCs)] exists and the solution of (3,6) is 

X(s) = [a(s) J"^B(s) . (3.7) 

Assuming m = 1 (extension of the following results to the multi-input 
case is straightforward), the transform of the jth state is 


and 



(3.8) 


(3,9) 


Cramer’s rule is used to evaluate x.(s), with D(s) ~Det[d(s)l and 
Nj(s) =Det[a(s)] with the column of 0(s) replaced by G, Due to 
the complexity of (3.7), it is no longer feasible to obtain analytic ex- 
pressions for the integrand but it may bo evaluated numerically. Since 
the elements of d(s) contain C(s), will have a branch cut along 

the negative real axis and the contour of integration may be deformed as 
shown in Fig, 1 1-6 giving 
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whei'o 2, ##•, N ai^c the poles of (3,8), and the residues are 

evaluated at the poles by 


D’(s^) D(s)-D(s^) D(s) 


with /\i - s - s^, (Since the poles are determined numerically by 

t 

iteration, D(s^) =0.) The poles due to the structural equations of 
motion may be assumed to be complex conjugates with s = a. + ib., 

X/ u i/ 

~ ~ Poles due to IKs) may be real or complex but for the 

following development, it is assumed there are N complex poles within 

the contour. Since Xj(t) inust be real, the integrand must be pure 

imaginary and therefore x (re^^) and x,(re are complex conjugate ex- 

J J 


expressions. With N.(s ) 

J i> 


N^» + iN^ 


N/2 

= t 


2e N . cos b.t-N.„sin 

L Ji i 


^00 

in b^t “ J Im|^Xj(re^^)J e ^^dr, 


( 3 . 11 ) 


S.(t) = X (t) + X (t) . 

■^r '^nr 


The incompress ible flow transient response of the three DOF section 

of Table II 1-3 will be calculated for a unit step input command to the 

flap 'li(s) = for U/b = 290 sec Figure III-l indicates 

that the bending mode has a subcritical damping ratio of Q = 0,03 at 

this airspeed. To study the effect of changes in airspeed on the non- 

rational portion of the response, ) ] , j = 1, 2, 3 is plotted 

*J — 1 

in Fig. I1I-3 for U/b = 200,290, and 350 sec . At time t. x ^ (t) 

’ -Jnr 

is given by the integral of the product of the function shown in the 
-rt 

figure and e , The value of (0) is proportional to the area 

-Jnr 

under the curves, and since all of the functions go to zero at r r= 0, 

lim X. (1) :: 0, In other words, the uonrational portion of the response 
—Jnr ^ 






INTEGRANDS OF NONRATIONAL PORTION OF RESPONSE 
DUE TO STEP COMMAND TO THE FLAP . M = 0 
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I 


(loi's iioL partioipaLi' in the motion charaottovl.stic o.(. an nnstaljlc i] u(. iorlnp. 
airi’oil, Not(f that all of tlio curves nhown in I'ii’,'. ili-h ar(> .•-aiioothly 
varying iuncttons o(' U/b, <!ven an tho alri’oil t’iuttors at U/h ~ ;u)0 


Tlio poles and residues required to evalute (a.li) were calculated, 
and the integral was evaluated numerically. The comyioiumt and total, trans- 
ient responses of the? plunge and torsion modes are shown in ill-d. 

In this case, the contoui’ integral about the infinitesimal circle at the 
origin in Fig. II-6 will be non/.ero, its value being the steady state 
value of X (t) due to the step change in f3 _ . These steady state values 

were determined from (3,6) by applying the final value theorem, lim sX (s) 

j 

lim x.(t) rather than by contour integration. The oscillatory harmonic 
mode superimposed on a^(t) and to a smaller extent on to 

the very lightly damped flap mode which is not shown. The nonrational 
portion of h(t) is 75 percent of the rational portion of h(t) at 
t = 0, while the corresponding percentage for a(t) is only 15 percent. 
As in the case of the transient loads, the nonrational portion of the 
response is characterized by n rapid initial decay followed by a slow 
asymptotic decay, the entire function being a monotonically decreasing 
function of t. Since the response of a mechanical system to a step input 
in torque must start at zero, the sum of the rational and nonrational 
portions should cancel the steady state value of , Hence, the small 
nonzero value of a(0) and the larger value of h(0) are attributed 
to numerical inaccuracies in evaluating the residues. 

The following comments are mado^wlth respect to Fig, 111-4, 


1, The oscillatory motions typifying flutter phenomena are due 

entirely to the rational portion of the response. If a method 
were available for modeling only this portion of the system, 
it would serve to describe the pertinent features of the 
flutter problem. This concept will be pursued in Ch, V. 


2, The effect of the nonrational response on the oscillatory total 
response would tend to complicate the determination of the damp- 
ing ratio of the rational portion. Techniques for determining 
damping ratios which do not address this fact may produce in- 
consistent damping estimates. This effect may be aggravated in 
cases with random structural excitations. It an estimate of the 
nonrational response were available, subtraction of this estimate 
from the total response may improve the damping estimates. 
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RATIONAL, NONRATIONAL, AND TOTAL RESPONSES DUE TO A STEP COMMAND 
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Cluiplor IV 


PAJJK APPUO;:i^MAOTS AND AUGMKYrj-D STATK MKTHODS 


To obtMin solution^: of Lh(? u(*ro(vLas t ic (M|u:\ I i ons of inoU.oii, K({, 

(2,1), it: Is lUHiessary to S|)ucify tho aorod ynamic loads, 1,, In Chajitro* 

II tho c'xacL analytic lr)ads wore i^ivon for two-dimoiis lonai i ncompross 1 bio 
and supcM’Sonlc flow, and i\ tochniquo for obtaining similar loads L'oi* other 
flow regimes was indicated. In Chapter III, these loads wore incorporates! 
ii\to the aeroelastic equations, and solutions were obtained by an iterative' 
search procedure. In the presoni; chapter, the use of Pad(5 approximants 
of the loaels to produce augmented state aeroelastic models will be stueiied, 
rile acivantage to be gained is the ability to j)crfe:)rm analysis with tho 
resulting <'onstnnt coefficient, ordinary differential eeiuatiems, Tlie 
jienalty paid to achieve this advantage, .is that higher order models must 
be manipulated. An implicit constraint of the technique is the minimiza- 
tion of the required number of augmented states to adequately represent 
the loads , 


INCOMPRESSIBLE TWO-DIMENSIONAL FLOW 

Augmented state methods for this flow regime derive from R,T, Jones ^ 
[291, [9], exponential approximation of Wagner’s indicial loading function. 

Many investigators have used the method, including Golami and TiUke | 30], 
Baird and Kelley [31], Dugundji [32], Richardson [271 , and Lyons et 
al,,'69l. Jones* approximation is 


kj(f) ?-■ 1 . 


(. 1 . 1 ) 


Garrick [70] noted that the linearity of the goveiniing equations allowed 
the calculation of arbitrary transient .lift, functions by the convolution 
i nt(?gra I 


Q(O)kjd) I- \ kj(t--t) 


( 1 . 2 ) 


I 


nr <‘qiu v;i lent ly , 


-iiU) ^ 

2n/^jU- 



l<|(t-T)Q(T)cli . 




SiiuM^ Tiiiplt^co I rimsf ormnt ion of (4,3) yields 


P(s) = 2rtPbu|o.^ -1- £[k^(t) jJq(s) . _ (4.4) 

.Jones’ approximation, (4,1), gives the transfer function relating Q to P 


as 


P(s) 


^-2 



2npbU 

0,5s -f 

0,2808s + 

0,01365 

Q(s) 

-2 

0,3455s + 



s + 

0,01365 


It is well known from linear. system theory that the functional relation 
given by (4,5) may be described in the time domain by the constant co-* 
efficient, linear, ordinary differential equations 


X 

-1 


-2 


P(t) 


= -0.01365A X, - 0.3455(-j) x„ + Q(t) 

b 1 b "2 

- 2jtDbu|o. 006825 x, + O.lOSOsA x + 0.5Q(t)L 

( - b ~l b “2 ; 


(4.6a) 

(4.6b) 


Garrick proposed the approximation k^(t*) - (t * + 2)/ (t ^+4 ) , but this 
function does not have a rational Laplace transform and the resulting 
ai^proximat ion to P(t) cannot be given by ordinary differential equations 
as in (4,6), Hence, in order to ensure the computational efficiency ob- 
tained by differential equations, it is customary to utilize approximations 
whose transforms are simple functions, 

— —10 

The rational transform in (4,5) may be evaluated at s - re and 
the resulting real and imaginary portions compared to the generalized 
riieodorsen function shown in Fig, 11-5, Figure IV— 1 indicates that the 
approximation is a good representation of C(s), especially in the right 
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hal r~ plane , The aKr^n'meuL de tuj’ioi'a tt*s a.s o 

rhis is primarily duo to the j)roxiinity of l;lu‘ 
s r: -o,()45r) and -0,d. 


I ncroas( 
i sol at: 0(1 


s beyoiul 
poles at 


T Luuro IV-2 is a comparison between the exact nonrationai 
etficieut shown in Fi^, 11-7, and the lift coefficient ‘<iven by 
for the motion ^iven by (2,49), The exact coefficient is j^iven 
last term of (2,50), while tlie approximate lift coefficient was 
by partial fraction expansion with Q(s) ^ (r/(s+a) +(o as 


lift c o- 
(4.5) 
by tJio 
obtained 


nr 




, -at ’ 

-23t/0,5e sinu)t* 


au) 


( 0 . 0455-a)^+cj 


r-0.0455t» -at»/ 0-0.455 . . \\ 

e -e ^cos cot’ + ^ sin cot’/J 


bto 


(0, 3-a) 


-0.3t’ -at’/ 0-0.3 -.A 

e -e (cos cot’ + — = — sin cot} 

\ CO 


■J! «. 


7) 


where a = lo = 0.1414, a = 0,0074999, and b - 0,10055, Tlie approximate 
lift matches the exact lift very closely even for this heavily damped 
airfoil motion with Q - 0,707. This would indicate that equations (4,5) 
and (4.6) may be used to calculate accurate circulatory loads for incom- 
pressible flow. 


A unique feature of the incompressible case is that all of the cir- 
culatory loads involve the single nonrationai function C(s), greatly 
simplifying the approximation problem. Equations (4,6) are in a form 
which is compatible with the structural equations, (2,1). The resulting 
model uses the augmented states 


and is givcui by 


X 

”P 







! 

j 




— (Ut— 




I 


I 0 0 


0 O x 


.0 0 I 



0 I 0 1 Tk 1 ro' 


(4,8) 


K J L >< 

1 2 ^ 


where 


M - riM 
s ' nc 


= + 0-5 


""s - (V + 0-5 

Tl(^) R 1^0. 006825 (^)^, 0.10805(-g)] 


-- 0 . 01365 (^)^ - 0 . 3455 (^) 

o b , 


For an airfoil with n degrees-of— freedom, . x is n-dimensional, while 
is two dimensional. The submatrices in (4.8) are dimensioned 
conformably with these vectors and the tctal dimension of the model is 
2n + 2. Since the ’inertia matrix’ of the left side of the equation is 
nonsingular, multiplication by its inverse gives the standard form used 
by control engineers 


X = PX I- G^u 


(4.9) 


-G8* 


I 



"X 


0 

I 

0 " 


• 0 ■ 

X - 

• 

X 

F = 

-(M* )’^K' 


(M* 

, 0, - 

(M» )”^Ci 


X 


- h 

^'2 

F 

P 


- 0 , 


Tlie elements of F are functions of (U/b) and the eigenvalues of 
F are approximate roots of the aeroelastic equations of motion. Figure 
IV-3 compares these eigenvalues to the exact roots of the section of 
Table III-3, From the close agreement between the exact roots and the 
approximate roots, it is concluded that the linear rational model of the 
incompressible two-dimensional section, (4.9), is interchangeable with the 
exact model, (3,4), for the purposes of engineering design. 

It is also possible to compare the frequency responses of the exact 

and approximate models. The frequency response of x^^ due to sinusoidal 

oscillation of u. is obtained from (3.4) by tabulating (x./u.)(ico) = 

i J 1 j 

Nj (iw)/D(iu)) . Similarly, the frequency response is calculated for the 

approximate model of (4,9) by tabulating the transfer function of 
(X./u.)(s) for s= ico. Figure IV-4 compares (h/p ) (iw) and (O'/p ) (iio) 
for the section of Table II 1-2. The good agreement between the frequency 
responses, especially in the range of flutter at o) = 70 rad/sec, indi- 
cates that the poles and zeroes- of the approximate model provide valid 
representations of the exact system. The dip in the amplitude seen on 
all of these frequency response plots in the range 70 < oj < 80 rad/sec 
indicates the presence of complex zeroes near the flutter poles. The 
location of these zeroes is critical to any active flutter control 
scheme and they will be studied in detail in the next chapter. 

To apply the above technique to other aerodynamic regimes, indie ial 
loading functions must be calculated so that the exponential approximations 
may be obtained. Although there is a significant literature concerning 
such functions [Lomax et. al, Ref. 10; Drischler, Ref. 711, their calcula- 
tion is laborious [e.g., Roddon and Stahl, Ref,.. 72], and the technique 
has not been used widely. 
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COMPARISON OF ROOTS OBTAINED USING THE GENERALIZED 
THEODORSEN FUNCTION AND R. T. JONES' APPROXIMATION 


TO C ( s ) AS A FUNCTION OF U/ba'. 


a 



.5 


Exad 

Approximate 




100 
0 

Phase deg -100 

-200 

-300 

(a) PLUNGE 

FIG. IV-4 COMPARISON OF FREQUENCY RESPONSES OF PLUNGE y'ND 
TORSION DUE TO FLAP DEFLECTION OBTAINED FROM THE 
EXACT MODEL AND R. T. JONES' APPROXIMATION FOR C(s ) 

-n- 

✓ 



50 100 

u, rad/sec 



i 


U, VK]^A^S PADK APPUQXIMAN^r MKTIIon 


The |)rol)lom.s associatcMi witli (:hr (tniculalion- f)f liullc. ial loading 
runetions led Vt^pa lo an alternative method [d9], [94], The avail- 

ability oi Uio aerodynamic loads ior simple liarmonic motions, from 
numerous well-developed techniques kernel function or finite 

element methods), suj^Kested the approximation of the loads by Pade approx- 
iaiants. A Pade approximant (PA) of a function is a ratio of two poly- 
nomials wliich approximates the function in some range of its argument, 
Baker [73] gives a thorough summary of the pro )*^rties of Pade approxi- 
mants for the case where a Taylor series expansion of the function is 
available. Tlie usefulness of PAs is due to the ease of the analysis of 
the resulting analytic, rational functions as opposed to the original 
function. If the original function is known only in tabulated form, as 
for the aerociynamic loads, a PA may still be obtained by fitting the 
rational fraction N(s)/D(s) to the tabulated values for s - Ik using 
(for example) a least square technique. 


The intent of the application of such approximants to unsteady aero- 
dynamic loads is to allow the aeroelastic equations to be solved for 
arbitrary motions, i,e., throughout the s-plane. This assumes the 
analyticity of the unsteady loads, a point of some confusion in the past 
due to the discussion of the generalised Theodorsen function. Neverthe- 
less, such approximat ions have been utilised, a primary example being 
the design of the active flutter control system described by Roger and 
Hodges [4], 


Vcpa suggested that the PA of the generalised force, in the 

i^^^ mode due to deflection in the mode could best be represented by 


Nil 

a s i- n S i- • • • I i\ s -i a 

q, .(s.M) = -2 1 ^ Nl-1 

N . , N-1 


Ms I- b,,s 


(• b. 


Ni-1 


(4,10) 


This is referred to as an |N,N+l] PA, 
Ihe correct steady-state value 


The PA is constructed to yield 
and may be modeled by N constant 


-73- 


t 


fmvL'ficiioiil , ordinary dif ;(:or(Mil; la 1 oquaLioiiH, 'Vhv hinii rroqut^ncy limit, 
may l)o choson to ^ivo tlu? plnton theory load ^ while in incompr(‘r.- 
sible flow with M 0, a^ yields tlu' none Ireiilatory , ^virtual maits’ Icnul, 
Vepa [dd] shows that the numerator eoef f icients , a., are determined 
by the denominator eoef f leicMits , and by eoastralnts (e.p;, , a ^ = 

dotermlnat ion of the PA for ^ involv(*s the deter- 
mination of the N constants, b^, i 2, N-i 1 , Hence, if ci^^(ik,M) 

is tabuliitod for more than N values of k, the b. may be determined 
by a linear least squares teclinique. 


Vepa comments that in order to obtain accurate l^As , a large number 
of reduced frequency values are required in the high fx’equency and the 
low frequency range. In his Ref, 33, tables of [2, 3] PAs arc given 
for plunge, rotation, and flap modes in two-dimensional flow for i\I r- 
0.3, 0.4, 0,5, 0.6, 0,7, 1.5, 1.75, and 2.0, Tiie location of the 

roots of the denominator of (4.10) are of concern since the PAs of un- 
steady loads are approximations to stable dynamic systems and should 
themselves be stable. Also, in the incompressible case, the loads are 
multiple-valued functions with branch cuts and, since the PAs are single- 
valued, they cannot yield valid approximations along the cut. Baker 
gives examples indicating that the poles of the PAs cluster along the 
cut and conjectures that as the order of the PA tends to infinity, the 
poles tend to a ^pole-distribution^ along the cut. Since the incom- 
pressible loads have a branch cut along the negative real axis, the poles 
of R.T, Jones^ approximation would thus be expected to be found there 
and this is indeed the case, the stable poles being located at s -0.0455 
and s : -0.3, 


This approximation technique would be an attractive method for the 

analysis of the aeroelastic system were it not for the high order of the 

resulting model. With each generalized load modeled by its own independent 

2 

[N, N+l 1 PA, the dimension of the model is 2n 4 Nn . For a three DOF 
section with |2, 3] PAs, the dimension is 24; (i structural states and 
18 PA states. On a realistic design problem, the dimension associated 
with the approximants may become ten times that of tlie original struc- 
Lurai model. The design obtained in Rt'f, 4 uscmI from IH t:o 27 structural 
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mockiji, whil(' tlie total (liiiion^ion was as hlfth as aoo. It would obviously 
bo dosirablo to ol)taln oompai-ablo accuracy with lower ordc'r ap|ir<»xlma- 
I U)us . 


C, ^r iK MATRIX 3^Ai)K APPU0XI MAN7 


Examination nf the 

transi(?nt 

responses obtained by inverse 

Laplace 

ti’ansformat ion of the PAs led Vepa 

to attcmi)t the aiiproxiimitlon o 

:f tJic* 

generiil i/.ed load matrix , 


by matrix I’As [73,], T'bc load 

matrix 

ve fates fo L as 

L 

2 2 ■ 
Pb U 

= 0(s,M) X(s) 

(4.11) 

and the matrix PA is 

o(s,m) = 

= P(s)R(i)"^ 

(4.12) 

where , 

Hi) = E 

i=o 

P s" 
i® ’ 

R(s) = ^ R.s 4- . 

i=o 



The minimum number of augmented states is given by N - 1, 

Q(i,M) = CPq + + P^s^][ls + R^]“^ . (4.13) 

The properties of Q(s,M) are dependent upon the eigenvalues of the 

matrix and the philosophy underlying the use of the matrix PA is that 

the individual loads q. may be approximated by suitable linear combina- 

^ J 

tions of ’sharecr eigenvalues. As in the scalar case, the numerator 

matrices are determined by and independent constraints, while 

is determined by a least squares procedure. The numerator matrices 

allow three constraints, one of which is the matching of the steady-state 

loads, Q(0 M), This determines P as 

o 
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Q (o,M)ir . 

SH O 


IMu* I’DiuainlnK two cons tj;a ini. ,s may bo used to oni’orco a match l)ctwo('n tlu* 
matrix PA :;acl the tabuiatcMl slmpla* harmonic loads at the anticliuitcd 
rcduc.iMl j;roqiuM)Cy o£ riutto]\ and with Q(ikj,, M) 




(■1, ir.) 




( 4 . 1 ( 3 ) 


Although the piston theory l:*mit is not enforced, the piston -theory loads 
are used with (4.16) to obtain a solution for K^, The simple har-- 
monic loads are given at m values of reduced f reqaciUl.y^^>!±eJLdlAig 




°2<h- 


- Qp,p.: ^ = 1 m • 


(•4.17) 


Equation (4.17) provides nm equations for the n unknown elements of 

R and with m > n a least squares solution is possible. The matrix PA 
o 

of (4.13) requires only one augmented state for each degree-of-f reedom of 

the structural model, and the resulting model has dimension 3n as 

2 

opposed to 2n •+• Nn*^ for the previous PA model. 

fhe incorporation of the matrix PA of (4.13) into the structural 
equations is facilitated by transformation to the state space model — 


X 1" X 1 Cl X 

-p p-p p- 


(4.18a) 


2 2 
Pb U 


X i- IPX r H.x 

"p 1“ 2“ 


(4.18b) 


win* !'(' 


p 


(i 

1> 


-L p - P.K I [p - \\\i 

a in P o o o 1 o 



1 


P 

a 


)V, K 
1. o 


1 P,^U 


2 

() 


"2 " - ''2'‘o> 


The matrix PA .nodol is [«iven by (2,1.) and (4,18) 



(4,19) 


The accuracy of the matrix PA will be illustrated in the remaindei' of 
this section. 


C , 1 Supersonic Matrix .Pade Approximants 

Simple harmonic oscillatory supersonic loads (2,67) were calculated 
for a three DOF section with the elastic axis at midchord and a 20 percent 
chord aileron (a == 0.0, and c = 0,6), The load matrix, Q(ik,M) was cal- 
culated for k = 0,05, 0,1, 0.2, 0,3, 1,0, 2,0, and 3,0 for the 
Mach numbers M =: 1,5, 2,0, and 2,5, The stead^^-state loads were cal- 
culated from the Ackeret— formula 


as 


4 


ViiM 


dz 

dx. 




0 4 -2(l-c)’ 


h 

0 i) (l-c^) 


X 

0 (1-c)^ -d-c)^ 

1 

‘ *, 





(4.20) 
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I’ho i)tsLoii tlu'ory loads aro dorlvcHl from tho startinn pressure Ap/ioU^ 









-4.0 

0,0 

“•0.16 
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l,iK 







0 

"1, 

0.1386 
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1 

O 

J-* 

a 

0.1 386 

"0,0486 


P 


L. 




— 


(4.21) 


1 lie matrix PAs were calculated foi’ M l.a, 2,0, .and 2,5 for the assumed 
flutter frequency - 0,2, and the resulting approximants are tabulated 
in Table IV-1, 

lor this case, the exact loads may be calculated for general values 
of s and compared with the equivalent loads calculated from the PAs in 
Table IV-1, Figure IV-5 gives this comparison for the loads c. . c 

^np = 2,0 for s = re^® with 60° 0 ;s; 150°. The two sets 

of loads are indistinguishable for jk| = r < 0.25 and generally agree 
to within 5 percent for r < 0.5, 

The PAs of Table IV-1 were used to calculate the eigenvalues of 
(4,19) for the section of Table III— 4, Figure IV— 6 compares those eigen- 
values to the exact roots of the aeroelastic equations from Fig. III-2. 
With b = 1.35m, and a^ = 333 m/sec, the reduced frequency is 0.11 < 
k < 0.23. While the PA was constrained to yield the correct oscillatory 
loads at k = 0.2, good agreement between the exact and approximate 
roots is seen throughout this range of reduced frequencies. 

From the discussion in Ch. Ill, Sect. A-2, it is known that the 
supersonic aerodynamics introduce an infinite sequence of poles of in- 
creasing modules to the aeroelastic system. Thus it may be anticipated 
that the eigenvalues of the PAs would provide estimates of these addi- 
tional poles. The eigenvalues introduced by the PAs are given in Table 

IV-2. Tho first column of the table gives the eigenvalues of R (scaled 

o 

by U/b), These poles are associated with the ’open— loop* aerodynamic 
medium as modeled by the PA,, The second column of the table gives the 
eigenvalues of (4.19) introduced by the PAs, and may be interpreted as 
the *closed-loop* poles resulting from the interaction of the structural 
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SUPERSONIC MATRIX PAD^ APPROXIMANTS FOR A THREE-DEGREE- OF-FREEDOM SECTION 
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FIG. rV-6 COMPARISON OF ROOTS OBTAINED USING GENERALIZED 
SUPERSONIC LOADS AND MATRIX PADE APPROXIMANTS 
AS A FUNCTION OF M 



dynamics with the norod^namic medium. 


Table IV-2 

KIGKNVALUbS 01’ SUPERSONIC PADE APPROXIMANTS 


M 

Eigenvalues of 

R (U b), rad sec 
o 

Eigenvalues of Eq, 
(4, 19) (rad/sec) 

1.5 

-690, -193 ± il47 

-844, 

-243 ± 1181 

2.0 

-1848, -484 f i331 

-1855, 

-480 ± 1333 

2.5 

-2784, -657 f 1448 

-2785, 

-655 ± 1407 


For this three DOF system^ the eigenvalues of the PAs are characterized 
by a complex conjugate pair with slightly greater than critical damping 
and a real root of larger modulus. All of the eigenvalues are stable 
and inci'ease in modulus with increas ing Mach number , This correlates 
with well-known results of piston theory, in which no augmented states 
are required to model unsteady aerodynamics at hypersonic velocities. 

Tlie eigenvalues are well above the bandwidth of the bend ing-tors ion 
airfoil section (o) < 100 rad/sec) and the slight difference between 
corresponding eigenvalues in the Table indicates that the modes of the 
PAs do not couple strongly with the structural modes. It is interesting 
to note that the complex pair at M r= 2,0 do not correlate well with 
the exact eigenvalue of the aeroelastic system at s - -1315 ± il501 
determined in Sect, III-A-2. 

C-2 Subsonic Matrix Pade Approximants 

Vejja has calculated matrix Pade approximants (PAs) for a three DOF 
section in subsonic flow and provided them to the author in a private 
communication [741, Tables IV-3 and IV-4 ])resent these PAs fer M = 0,3, 

0,4, 0,5, 0,6, and 0,7, Tabulated oscillatory aerodynamics at k 0, 

0,3, 0,4, 0,5, 0,6, and 5,0 were used to construct the ai)proxlmants 
ot ’rabie IV-3 which assumed k^, -- 0,5 while the approximants of Table 
IV-4 include the reduct?d freciuency k 0,01 and assume k^. - 0,4. 
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Table IV-3 

SUBSONIC mXRIX FADE APPROXIIViANTS FOR A THREE-DEGREE-OF-FREEDOM SECTION (K = 0.4) 
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Table IV-4 






I 


Vcpa'.« sign conventions differ from those used herein In thf 4 t: h, P are 

positive upwards, and P,M^’ nrc' positive for trailing*-edge upwards rotation of 
th(' aileron. The pitching moment is calculated about midchord and a quarl(>r- 
chord aileron is assumed. The oscillatory aerodynamics were calculaled by 
Vepa using a kernel function program described in Ref, 33. 

While no calculations are available to check the validity of these 

PAs for arbitrary complex values of. s, they may be compared to known 

simple harmonic oscillatory loads. Figure IV-7 compares c., , c,„ . and 
_ Zh '"a 

Cj^p for s = ik with the Tables [75], [76] based on Timman and Van 
deVooren’s Mathieu function solution [12]. The PA of Table IV-3 gives 
an excellent match ejccept for ReCCj^^) which is 10 percent low at k = 0.3, 

The approximants of Table IV-4, with the low frequency point k = 0.01 
included, show a deterioration in agreement with the accepted values. 

The constraint imposed, on the latter PA by requiring agreement at k = 0.4 
is clearly seen, however. It is concluded that with appropriate checks, 
the matrix PA technique can provide a good augmented state model of 
oscillatory subsonic loads which are also valid approximations in the 
vicinity of the ito axis. 

The eigenvalues of the R matrices of Tables IV-3 and IV-4 are 

o 

given in Table IV-.5, All of the eigenvalues are distributed along the 
negative real axis except for the PA of Table IV-4 at M = 0.3 and 0.4. 
Evidently, the inclusion of loads at k = 0.01 causes the approximant 
to develop complex conjugate roots with the resulting deterioration in 
agreement at higher values of k shown in Fig. IV-7, Henceforth, only 
the approximant of Table IV-3 will be used. 
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FIG. r\"-7 CO.MPARISON OF SUBSONIC OSCILLATORY LOADS CALCULATED USING MATRIX FADE 
APPROXIMANTS AND THE MATHIEU FUNCTION SOLUTION AT M = 0.5 



1 25 Mathieu function solution 

i o Tabie IV-3 (k, = 0.5! 



Tab.l(‘ IV-fi 

KKiKNVALUEi!) OF \\ FOU HU13SONJ:(’ MATItlX ]>AUK Al^J^KOXIMANTS 


Table.' IV-:3 


Table IV-/} 


0,09504 

-2,807 

-5.480 

o.ioaf) 

-2.560 

-3.250 

0, 1218 

-1.86:5 

-:5.0H4 

0, 1091 

-1.244 

-1,798 

0.087:M 

-0,5887 

-1,137 

0.f)2374 

-3,541 ± 

i 1.560 

0.02262 

-2,533 ± 

i 0.6013 

0,02212 

-1.263 

-2, 34 3 

0,02105 

-0.6990 

-1,782 

0,04512 

-0,3401 

-1.221 


£. STATIC DIVERGENCE 

If the incremental moment generated by airfoil pitching is greater 
than the restoring moment of the torsional spring the airfoil is 

said to be statically divergent. The divergence velocity is given by 
[52, p, 193] 


Static divci’getice in incompn 
DOF pitch equation 


2 

2b\h - a) ^ 


■5sible flow may be studied using the single 


1 ^ s u'( s) iK s) 2TfPb 'u( y+a)c( s)[ b(|-a)s i U] >:( s) 

Substituting (4,22) and assuming 5c /Jxi 2Tf 









1 


) 






b(| !-;))« t U 


' '‘‘cj 




(4.24) 


Since static divergence is a low frequency phenomenon, the s and s terms 
may be neglected, giving 


a(s) = 0 . 


(4.25) 


A pole of the aeroelastic system occurs at values of s for which the 
coefficient in (4.25) is zero. Since .C(s) is purely real only on the 
positive real axis, poles can only occur there. Also, along the positive 
real axis, C(s) decreases monotonically from a value of 1.0 at r 0 
to 0.5 at r = Hence, a pole cannot occur for U < and for U > 
only one real pole can occur. This mode produces the motion of the diverg- 
ing airfoil and occurs in addition to the 2n structural poles. 

The occurrence of this divergence mode may be studied by locating 
the poles of the system in the s-plane. The exact system model of (3.4) 
or the Fade model of (4.9) may be used to locate these poles. The airfoil 
described by Table III-3 was modified to yield a divergence speed close 
to the flutter speed by setting g, = 5 and a = 0. Table IV-6 compares 
the pertinent roots of the exact and Fade models as a function of airspeed. 
The plunge and torsion mode poles are given for both models and the aug- 
mented state pole with largest magnitude is given for the Fade model. The 
divergence speed of this three DOF section at 1.13 is indicated 

for the exact model by the emergence of an additional real pole on the 
positive real axis. Since the K;T. Jones approximation to C(s) is 
accurate at low i rcquenci(?s , the Fade model gives a valid approximation to 
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Stahl [72] using strip theory and uuKinontud states. They also found the 
static^ divergence mode ro be given by the* augmented state with largest 
magnitude. 


Table IV-6 

STATIC DIVERGENCE IN INCOMPRESSIBLE FLOW 
(poles in rad /sec) 


u 

Exact Model 

Padd Model 

1.10 

7.407 ± i 73.51 
-36.10 ± i 48.14 

6.974 ± i 73,29 
-25.61 ± i 31.79 

-0.6253 

1.11 

7.825 ± i 73.18 
-37.02 ± i 48.04 

7.314 ± i 72.96 
-25,85 ± i 31,80 

-0.347 

1.12 

8.150 ± i 72.86 
-37.94 ± i 47.95 

7.650 + i 72.63 
-26.09 ± i 31.81 

-0.0652 

1.13 

8.470 ± i 72.55 
-38.87 ± i 47.85 

+ 0.1885 

7.982 ± i 72.31 
-26.33 ± i 31.81 

+0*2431 

1.14 

8.786 + 1 72.23 
-39.80 ± i 47.75 

+0.4853 

8.309 + i 71.99 
-26.58 ± i 31.82 

+0,5702 

1.15 

' 1 

9.098 ± 1 71,92 
-40.73 ± i 47.60 

+0.8172 

8.632 ± i 71.68 
-26.82 ± i 31.83 

+0.9156 
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ACTIVK C:o>miOL OF AKUOKI,/\STIC^ SYSTFMS 


Tli(» aflvanc(‘s mado diuvlu^ Uk^ past doead(' in the rt» I i abl 1 i. L y and 
acceptability rsf active control tochniquos as a})pliod to the stability, 
control, and navigation functions of aircraft have induced similar 
advances in the aoroelastic design of aircraft, (iarrick [77] provides a 
synopsis of this activity, while Table V-1 lists tlu? categories commonly 
ascribed to chls control configured vehicle (CCV) technique. 


Table V-1 

CONTROL CQ1!FIGURED VEHICLES DESIGN CATEGORIES 

augmented stability (AS) 
maneuver load control (MLC) 
ride control (RC) 
fatigue reduction (FR) 
gust alleviation (GA) 
flutter mode control (FMC ) 


The different categories in the Table have traditionally been characterized 
as affecting either the low-frequency, rigid body response or the high 
frequency, elastic mode response. This ^bandwidth separation^ in the CCV 
functions has become a moot subject as larger and/or more flexible air- 
craft are designed and the analysis of the interaction between the rigid 
body and elastic modes is becoming commonplace. 

To control the aeroelastic system, it is necessary to apply a control 
force or torque. Although nonaerodynamic controls have been considered 
(78, Buchekj, current designs use aerodynamic surfaces to produce the 
control leads, The B-52 Load Alleviation and Mode Stabilization program 
(LAMS) (1| utilized accelerometer measurements at the location of the 
control surface to produce augmented damping of subcritical structural 
response modes, .Tn the subsequent I3-r32 CCV program [ r9 I , [2], all of 
the items of Table V-i we>rc incorpoini ted . Roger (uid IIt)dges fd] document 



iho llutuor wuu\v control .system tests of tliis program, wjiich is the 

only flifvht Icstocl FMC system in exist once, TIris FMC ciesin’ii [2] utilised 
auniuent(Hl state l\ulo appi’oxxinauts for tlie loads. Grosser et al, [3], 
outline the C-5A active lift distribution control sy.stem wliicli incorporates 
MCL, FH, and GA systems. Wind tunnel studies of FMC systems are reported 
by Sandford ot al, [3], and Ilaidl et al, [80], The former test was de- 
sijjjned usinj^ Nissim^s aerodynamic ener^^y technique [3^] and encountered 
difficulty in stabilising a leading edge control surface. The latter 
test studied wing-store flutter and incorporated a FMC system designed 
using oscillatory loads and a simple damping control law. 

It is obvious. that an analysis technique capable of treating FMC 
can also be applied to the less demanding tasks of AS, ,RC, FR, and GA, 

Also, the FMC problem provides a definite design goal-stability , whereas 
the other CCV categories have more subjective design criteria. Thus the 
FMC problem has been the subject of a number of analytical studies,. . 

Turner [38] used a modified j>-k method with oscillatory loads to obtain 
a model amenable to modern control techniques while Dressier [39] 
used a series expansion in s for the loads to obtain an. augmented state 
model and applied modern optima'* control methods. 

This brief review of the literature of FMC focuses attention on the 
key role played by the choice of the aerodynamic model. Most of the 
above studies v/cre conducted by obtaining a model described by linear, 
constant-coefficient, ordinary differential equations. The complexity 
of the various aerodynamic modeling techniques ranges from no augmented 
states [38] to well over 100 augmented states [4], It is significant 
that the only flight tested system, the B-52 CCV FMC system, used the most 
complex aerodynamic model. Garrick [77] compares the predicted flutter 
characteristics of the analytical model, the wind-tunnel model, and the 
flight tes t results of this program. The general trends of the damping 
of the flutter mode were predicted accurately, although the predicted 
flutter speed was off by 10 percent. Thus, there is room for improvement 
in the modeling of aoroelastic systems. Desirable characteristics of 
i mp r o ved mod els include : 


( 


(1) imprcwcd accuracy in prodictiuK arbitrary transient response, 

(2) minimi /-at ion ox' elimination of auM;inented states required to 
model the loads, 

(d) maintaining’ the simplicity of ordinary diff erent iai equations 
tor the mociol , 

(4) applicability of modern control techniques to system synthesis, 

(5) applicability to flight test results from the vehicle for 
which the active control system is to be designed. 

The last item is stressed since it implies the possibility of tailoring a 
system to a vehicle during a flight test program. It., might be hoped that 
future FMC systems will not require the degree of analytical study of the 
vehicle which was available to the B-52 CCV pi^ogram. The design technique 
which will be developed in this chapter addresses itself to the above items, 


A. CONTROL OF DISTRIBUTED PARAMETER SYSTEMS 


The structural elements comprising an aircraft are three-dimensional 
elements (wings, fuselage, empannage, tail) whose dynamic behavior is des- 
cribed by partial differential equations with appropriate boundary condi- 
tions, Bisplinghoff and Ashley [52] formulate the aeroelastic equations 
in terms of operators as 


(S - “o - " % 


( 5 . 1 ) 


Where J ai’e structural, aerodynamic, and inertial operators, 

q is a generalized displacement, and is a known disturbance force. 

Depending on the formulation adopted, (?^, or J may be algebraic, 

differential, or integral operators. For instance, the structural operator 

2 2 / 2,2 

for the bending displacement of a simple beam is = 5 (Flo /3y )/dy . 
Neglecting shear deformation and rotary iziertia, the uniform slender-beam 
differential equation is 


PT 1. 

dy di 




(5.2) 
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tho S{)lullnn of (T),2) 


For cantilever bounciai'y (anulitions and 0, 

can be expi’ossecl as [7| 

CO 

w(y.t) = £ <t),(y)§.(t) (9aa 

k=i *' 



(t) =: B, sin oj, t + C, cos lo, t 

k ^ ^ k k k k 

2 El 

a = — . 

m 

The natural frequencies u) are given by the solutions of the equation 

K 

1 + cos^~ ^ cosh X> - 0 , 

This example illustrates the key concepts embodied in the study 
of the control of distributed parameter systems. This field embraces 
the study of linear operators defined on a Hilbert space and seeks solu- 
tions to ’optimal control' problems specified by an appropriate per 

formance index. The distinguishing feature of such problems is the 
infinite dimensionality o.. the solutions- (or the elements of the space). 
This effect is clearly evident in (5.3) where the solution is described 
by an infinite sequence of normal modes. Much effort may be expended 
in establishing the existence of bounded inverses of the operators since, 
in this event, the solution may be uniformly approximated by a finite 
sum of 'normal modes’. (A given function is uniformly approximated by 
a sequence of functions if the approximation invariably becomes better 
as additional elements of the sequence are incorijorated . ) This is the 
basis of the well-known method of truncated normal modes in structural 
dynamics problems. Since the infinite sequence of orthogonal 'In vacuo' 
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mo(lc\s oJ: (5.3) span Uho solution sjuico^ tlu* solution to lIk' .'U't.uai 
problcMii, with / 0, may bo ol)tnino(l l^y pro.jr'tr t i ni;- th(' appliful 

clist r ibut(Mi i’orc(' onto tlu\scj M^asis vectors^ 'ilio rt.-sultluK 

^CMU' !*al i/.o(l to t’cr os , a I’o fi;i vc*n by 

o 

In aeroelastic problems, the applied force is composed of forces due to 
motion of the structure and 

where is the generalized force in the ith mode due to deflection 

of the structure in the mode. 

In the above example, the infinite dimensionality of the solution 
is explicitly indicated by the partial differential equation, (5,3), 
describing the structure. ..The examples treated in the previous chapters 
involve pitch and plunge of two-dimensional typical sections which may 
be regarded as representing the first bending and torsion inodes of a 
three-dimensional flexible wing. The dynamics of such t 5 q)ical sections 
are described by ordinary differential equations. However, even, these 
cases require infinite d iinens ional solution spaces since the aj^pllcd 
loads are themselves solutions of partial differential equations (e.g., 

Eq, 2,6). In the former case (elastic structures ) the spectrum of the 
structural operator contains an infinite sequence of discrete eigenvalues, 
while in the latter case (typical section) the aerodynamic opcu’ator 
may have a continuous spectrum, as typified by the branch cut of C(s), 
or it may be discrete, as In the case of two-dimefhsiona 1 supersonic 
loads, Note that the elastic structure problem Involves the solution 
of two distributed parameter systems and the solution sp(.‘C‘lrum will he 
doubly inl'inile, A 1 Uiougli tiu’ matlu’mat icnl desc r* i p (. ion of’ Uu* str’uctur<‘ 
and tlie aerodynamic medium thus apjiear to bt* oti an tnjual ionting, 
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the HLructur<' wh uMi i .s iiivnriiib^y viewod ns Dip Lo bp pon LiV) 1 1 pd , 

WnuK 'runj^ [^2] providod a framowoi’k f(^r Ihr. sLudy of disLribuifd 

pai’nmpLor control, and ouL 1 i nod , Xho onte^orios of (1) d i s tr i butf'd 
input control, (2) boundary input control, and (:i) total Input control. 
They oxtondod the conco|)ts of controllability and observability, which 
were developed with regard to finite state space, to the infinite 
dimensional case and examined the problem of existence of solutions. 

Taking the view that an ’optimal control’ should be defined with respect 
to the complete solution of the mathematical problem, they were Led to 
performance indices defined on direct sums of Hilbert spaces. Problems 
posed in this vein have proven unwieldy with the examples considered 
usually having only one spatial dimension. References [43] through [46] 
illustrate the theory applied to the one-dimensional heat equations, 
while [47] and [48] study the one-dimensional hyperbolic equation and 
wave equation. It is of interest that several recent references [46], 
[50], and [49] address the more modest goal of ’stabilisation’ of dis- 
tributed y:>aranieter systems rather than seeking an ’optimal control’ in 



In assessing the relevance of distributed parameter control theory 
to the aeroelastic problem, it must be noted that none of the three 
categories of control given in the last paragraph correctly describe 
the problem. The control force available in the aeroelastic problem 
is the pressure distribution caused by control surface deflection. 

It cannot be considered a distributed input since it is a one-dimensional 
function of the control deflection. NeLther is it a boundary control 
for the structural partial differential equation since the boundary 
condition associated with the surface deflection relates to the aero- 
dynamic equation. Hence a broader formulation is to address the prob- 
lem properly. 

Jian and Ching Yuan [49] have presented such a formulation. They 
model a distribuleci parametc*r system with an ordinary feedback controller 
as 


- 98 - 


I 


m( y) 


p\ ^ pi 

— i ^ — Cw f Bw i* Aw 

Bt Bt 


-Gx 




dx 

dt 


JX h 


S^v S ^ 


(5.4b) 


2 

where A, B, and C are matrix operators defined on a Hilbert space L , 

and G is a bounded operator mapping the n-dimensional vector space 
2 

into L , G is thus the operator relating control surface deflection, 

X, to a pressure distribution over the surface. The ordinary feedback 

1 ? 

control is derived from the n-dimensional vector x. S and S are 

2 n 

’observer-operators* mapping L into R , In [49], the model of (5,4) 
is analyzed from a rigorous Hilbert space standpoint. The stability 
of the system with feedback control is studied and several perturbation 
theorems regarding the spectrum of eigenvalues are proven. Finally, 
the validity of truncated normal mode approximate solutions is verified. 
Unfortunately, no examples are given in [42], 

The use of the truncated normal mode method lends insight to the 
concepts of controllability and observability of distributed parameter 
systems. In the context of aeroelastic wings whose motions are measured 
by * point sensors * (e , g, , rate gyroscopes , accelerometers’, etc . ) , an 
aeroelastic mode will be unobservable to a sensor placed at a node of the 
mode (i.e,, if the measurement distribution vector Is orthogonal to the 
modal eigenvector ) , *S imilarly , an aeroelastic mode is uncontrollable 

by an aerodynamic control surface if the generalized aerodynamic force 
in the mode due to control deflection is balanced by the remaining 
elements of the aeroelastic equation, (5,1), In control theoretic terms 
this implies that the control distributor vector is orthogonal to the 
reciprocal eigenvector of the mode. The rigid two-dimensional sections 
unalyzt3d herein are certainly observable if both displacement and angle 
sensors are employed. The controllability of such sections with respect 
to leading- and/or tralling-edgc controls will be examined in the next 
section. 
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! 
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n-, CQOTllQUABXUTY ANB QB>SKUVABILITY OF AKUQEMSTIC MQDKS 

Controllability of tluj linear, constant coefficient, finite dimen*- 
Hionnl sy.stcm 


X = FX i- G^u 


(5.5a) 


y - HX 


(5.5b) 


was examined by Gilbert [82]. The dimensions are 

X " N-dimensional state vector 
u = m-dimensional input vector 
y = p-dimensional output vector 
F = NxN matrix 
= NXin matrix 

H - pXN matrix. 

If F has distinct eigenvalues, the transformation X = Zz where the-. 

columns of Z are the eigenvectors of F , gives 


z = Az + (Z ^G^)u 


y = (Hz)z . 


(5.6a) 

(5.6b) 


The elements of the diagonal A matrix are the system eigenvalues, 


i = 1, N. Laplace transforming (5.6) gives 

\ 

y(s) = y(s)tl(s) 


with 


3l(s) = (HZ) (si - A)~^(z"Vt^) 


N A . 


(5.7) 


(5.8) 


where is a pxni matrix given by the vector outer product 




(5.9) 




(ir/) . is thv i column o C HZ and (Z G.) is l\u^ i iii rf>\v (»1 Z '(; . 

» 1 1 i . r 

'riie modo r(‘))rcs(*nLtul by will be* uncont ml 1 al) I c U ihc i ,m\v ol 

-1 ^ 

Z Gj is zero and 11 will bi* uno])servabl e 11 Mic i i li t'oluinn oi MZ u. 
z('i‘(), K(|Ualioii (b.f)) in^JU^a,U^s that in eiliior of t.iK-ie case's, A. o, 

aiul the sysli'iii Iransfor runcM. ions }*iv(‘n by l((s) will nol (ton la in I he 
polo at Xp in other weirds, l.f X. either unin)nl ro I Lable nv un- 

observable^ llien evoz’y liainsl’er 1‘unellon in Jj(s) wi ll have* a /.cm al 
X^ and a pole-zero cancellation will occur. The fact tliat a pole- 
zero cancellation has occurred is not sufficient information to deter- 
mine whether the system is uncoutrollable or unobservable, Tliis must 
be determined by examining' the rows and columns of the input and output 
matrices. These observations regarding' the relationship of pole-zero 
canv'cllations and controllability and observability are the basis of 
the design technique used for the B-52 CCV flutter i ule control system 
[2], Control surface positions and sensor locations were chosen to 
achieve the largest separation between the flutter mode and the nearest 
zero. In a realistic design situation, exact pole-zero cancellation 
may not occur but a near pole-zero cancellation inay indicate that the 
required control power will be excessive. Also, near cancoliation fre- 
quently leads to severe sensitivity problems. 


£. CONTROLLABILITY OF A TWO-DIMENSIONAL TYPICAL SECTION 

Since the aeroelastic mode shapes of flexible wings vary contin- 
uously as functions of velocity and dynamic pressure, it may be anti- 
cipated that observability and controllability problcnns will occur at 
discrete values of these parameters, if tlicy occur at all., Tlie typical 
sections analyzed in the previous chapters are observable if measurements 
of h, a, and y are assumed. Tims, the cont rol labi li I y of the sec- 

tions may be studied by examining the transfer functions of the Pade 
ap])roximant augmented state models (MLade models*) given by (4,9) for 
the incomj)r(ms ible case or ijy (4,19) .Cor the compress i.b le case, 

'Cable V-2 gives the parameters defining a nominal case Cor tlu' 


four DOF section of Fig, 11-1 in incompr(*ssii)l(.^ flow. 


Ti>V)it) v-a 


N0MINA1.I PAIIAMKTKUS FOR A FOUR“OKCiRKK-OF-l''RKKI)OM 
SICCTION IN INCOMPRKSSIHl.F Fl.OW 


''a 

100 rad/ sec 

c 

0.0 

‘•'h 

r 50 rad/soc 

X 

a 

:: 0.2 


500 rad/soc 

2 

= 0.25 

(i) 

r 

500 rad/soc 

"p 

z: X r. 0.0125 

r 


t.* 40 

2 

-r r"*^ - 0.00625 

r 

a 

= -0.4 

% 

r 


The leading- and trailing-edge control surfaces span 20 percent of the 
chord and have natural frequencies five times the torsion mode frequency. 
A small viscous damping has been assumed to stabilis'.e the flap modes. 

The remaining parameters are identical to those of lable III— 3, 


It is well known [IS] that the frequency ratio has a strong 


influence on flutter characteristics. Therefore calculations made with 
the Fade' model (4.9) for 0, 0.25, 0.50, and 0.75 are pre- 

sented in Fig. V-1.. A subprogram was written which iterated to determine 


the value of ^/bw^ at which flutter occurred. The transfer functions 


of (h/p.)(.s), (a/P.)(s), (h/r^)(s), and (a/r^)(s) were determined at 

c c 


this value of U/bu)^ and at _+ 25 percent of this value by the method 


of [831. Figure V-1 indicates the. variation of the poles and /.erocs 


of these transfer functions as a function of U 


/'.) increases from 0 to 0.75, the value of U/l)..' at flutter decreases 
IF a ^ 


from 3.41 to 2.13. The variation of the /.eroes is of special interest 
since they determine the controllability characteristics of the section. 
It is evident that leading edge flap control will not experience any 
controllability |iroblems, since the locus of /.eroes of lioth ()g ) (s ) 




Ai\D (c. /o? IiN INCOMPRESSIBLE FLOW 



S inco 


and (O/V )(s) toad to rail oulsido the range u)^ < < (,^ , 

' ui h a* 

classical flutter of a typical section occurs will) a merging of fre- 
quencies within this range ^ a pole-zero cancellation evidently will not. 
occur- for .loadlngredge control of a two-dimensional section . 

The situation for the trailing-odge control surface is much differ- 
ent , with the locus of zeroes of (h/f^^)(s) and (a/p )(s) tending to bo 
in the range < to < to^. In fact, at ”• 0.5, there is a near 

pole-zero cancellation in both degrees-of-f reedom at U/bu> - 2,84, 

Ifence the flutter mode is nearly uncontrollable at the flutter speed 
for this section. This explains the choice of to^^ -- 50 rad/sec for the 
nominal case in Table V-1, It represents a ’worst-case” design problem 
and parameter variation studies about this configuration are of interest. 
The physical cause of the uncontrollable mode can be explained by noting 
that for this specific set of values of the parameters of Table V-1 
and the structural and inertial forces and moments on the main 

section cancel the incremental lift and pitching moment due to flap 
motions when the section is oscillating in this * uncontrollable' 
mode. 



Figure V-2 gives the modal composition of the eigenvectors of the 
nominal section (Fig. V-lc) at U/bco^ = 2,13, 2,84, and 3,55, The compo- 
nents of the eigenvectors are presented in complex phasor form and are 
referenced to the plunge mode, h, which is normalized to unity. The 
uncoupled flap modes (<i = 500 rad /sec, ^ = 0,1) at s ^ -50 _+ i 497 
rad/sec have been modified by the coupling, giving a higher frequency 
mode at a; ^ 590 rad/sec and a lower frequency mode at co = 265 rad/sec 
Tile higher frequency mode is predominantly a trai ling-edge flap mode 
while in the lower frequency mode the leading-edge flap predominates. 
Also, a significant reduction in the leading-edge flap mode damping has 
occurred which explains the necessity of the viscous damping Q . The 
remaining two modes are the (highly-coupled) bending-torsion modes, 
otie of wliich becomes the flutter mode. Additional studies of two DOF 
bending- tors ion sections (not shown) indicated a minimal Influence of 
llie flap dynamic coupling on the characteristics o f— t lie flutter mode. 
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' ' - U/bu 
Mode' 

2.13 

2.84 

3.55 


P 

P 

< 

P 

ff 

Trailing- 
edge flap 

^ V" ■ 

a 

y 

a 

a 


s = -46±i584 

s = -60±i592 

s = -74.6±i603 

Leading- 
edge flap 

a 

'' p^ 

a 

a 

^ p 


s = -6.8±i272 

s = -11.2±i265 

s = -17. 8± 1254 

Torsion 

/ 

a 

s = -3.7±i96 

y^p\. 

a 

s = -0.4±i76.3 
(flutter) 

a 

s = 16.9±i67.9 

Plunge 

a 

h 

y 

a 

y 

a 

y 


s = -4.7±i54 

s = -12.4±i66.4 

s = -36.2±i61.8 


FIG. V-2 MODAL COMPOSITION OF THE EIGENVliCTOHS OF A 

FOUR DKGREE-OF-FHKKDOM SECTION VERSUS 

U/bo', IN INCOMPRl-:SSIULE FLOW 
a 
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This is due to the very low inoriias of the flaps rolativo t.o idu* main 
section . 


Figure examines the controllability of t!io nominal seetir)n 

2 

with respect to variations of the parameters: g, , a, r , r Cr ), , 

(./ i f 

)» Co(C )» )• ratio, p, is the only quantity 

fi Y f y i' T 

which is directly related to flight condition. The remaining parameters 
are related to structural and geometrical properties of the section. 
Variation of g, and the parameters related to the main section 
r^, a) have a strong influence on the controllabilo^t^ of the section. 

As- might be expected, variation of the parameters lolated to the flaps 
have a small perturbing effect on the controllability. The behavior 
of the zeroes associated with the leading-edge flap remains unchanged 
for all of the variations of Fig, V-3, Figure V-4 shows the effect on 
the critical flutter mode at flutter due to variation of the trailing- 
edge flap chord, c. Variations in c *detune* the uncontrollability 
condition. 


Thus, from considerations of controllability, the leading-edge 
flap has advantages over the trailing-edge flap for active aeroelastic 
control purposes , This advantage is offset by (1 ) the large destabili- 
zing hinge moments which the leading-edge flap must carry; (2) the 
associated power required to move the f lap , and (3) the violation of 
the aerodynamic shape of the lifting surface in the critical leading- 
edge area. In addition, proper design of a trailing-edge controller 
may achieve the objectives without encountering a controllability 
problem. The successful flight test of the B-52 CCV FMC system indicates 
that this is possible. 

It should be noted that the desirability of leading-edge control 
is not so obvious on finite wings which have a sequence of structural 
modes, w . Figures V-1 and V-3 shows that the leading-edge control 

K 

does not encounter controllability problems because the zeroes of the 

relevant transfer functions remain outside the range (!■ < a' < u.' . 

h r/ 

This reasoning fails when applied to finite wings since then the zeroes 
may cause controllability problems with other modes. 
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FIG. V-4 


HF FFFECT OF THE TRAILING-EDGE FLAP CHORD, c, 
XT mutj c'ttit'T’pr MnnE AND ASSOCIATED ZEROES 
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Controllability of the throe DOF section in (•nmprossiblo flow may 
be studied using the Pade model of (4.19), Figure V-h indicates that 
the subsonic section of Table V-3 has much the same controllability 
problem us the incompressible section. 


Table V--3 

THREE DEGREE-OF-FREEDOM SECTION PARAMETERS FOR SUBSONIC -FLOW 


~ rad/sec 

= 317 rad/sec 

P 

b 4 ft 

a 1000 ft/sec 
u - 40 
c r 0.5 
a - -0.4 




0.2 

0.25 

0,0125 

0.00625 



Figure V-6 shows the locus of zeroes of the (Q'/B )(s) transfer 

c 

function of the three DOF section of Table III-4 in supersonic flow. 
The lack of sensitivity of these zeroes to Mach number, and the fact 
that the (h/p^)(s) transfer function has no complex zeroes near the 
flutter mode indicate that controllability of two-dimensional sections 
in supersonic flows is not a problem. 
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POLES AND ZEROES OF A THREE DEGREE-OF-FREEDOM SECTION IN 
SUPERSONIC FLOW CALCULATED USING MATRIX FADE APPROXIMANTS 



. 1 ^ 


n . AKu oKiws^ric controi, basku on 'nil'; co^NO’i n* o k 

AKiiODYNAM K! KNKIU IY 

N is si. Ill 1 Hi. I (U'vol()]u?d a t(‘aliniqut? ol aaiivo riul!:ei’ su))j))H:»ss Ion wiilcl 
is ba.Hod upon c.onsldtu'atiojis oi tlio (uioi^vy ri‘qui.r(u) in sustain slinplf* 
harmonic osci’lalions o.i a two DOK typical section. If the sir;n of tins 
cnerny is positive', tndicatiUK t'Uorp;y must l)o supplied to tlu^ section to 
maintain the oscillation, the section is stable. A negative onerpy would 
indicate that the airstrcaiJi was supplying ener^.^y to the section and it is 
assumed that the section would flutter if released. For forced simpler 
harmonic oscillations of the section the aerociynamic energy is riven l>y 






(•^>. 10 ) 


n n ; 


The complex vectors 4 i^^ are generalized modal coordinates asso- 
ciated with the aerodynamic energy, and the elements of the diagonal 


matrix 


are the real eigenvalues of the Hermitian matrix 


C-(Q2 -I- Qo) + 


The matrices Q and Q. arc the real and imaginary parts of the simple 

X ^ 

harmonic aerodynamic loads. P, being a quadratic form, will be positive 
definite if > 0, i - 1 , . , . , n and thus the section will be stable, 

Nissim noted that this ’stability criterion’ was dependent only U)K)ii the 
aerodynamic loads 4 iQ^ and did not involve the structural parameters 

of the section (^t , etc.), Ap])arently the stability of the section 

could be delei'inim.'d without regard to tlie structural dynamics of tdic sec- 
tion. Nisslm argued that this was a desirable iorinulatiou because of the 
wide variati()iis in flight conditions which an aircraft may experience. 

In order to achieve active control of a lluttering sect, iini, Nissim 
postulatc'd Mie c^ontrol. law 


4 


(r).xi ) 








With tho Qa]) doIMoc* ti(>ns uxprcHSCKl as functions of h .'incl ry, the loads 

duo to flap dof IcHilions could be calculator} anrl atldod to and 

The stability of tho section could then be dotenninod by examining tho 

sign of X and \ . This stability ci^iterion had to bo checked ()vox" a 
X 

range of reduced frequencies, since the reduced frequency of flutter is 
not determined, Nissim determined the ’optimized^ values of C and C 

X /i 

for a section with leading-^ and trailing-odgo controls as 


C 


1 


-0,05 -1,7 

0.5 1.0 


0.45 0.2 

# 

-0.5 1.0 


Nissim [81] also studied control with only a trai ling— edge flap and found 
that it was barely possible to ensure the positive clef initeness of [^As] 
over a range of k. Furthermore, the design was sensitive to variations 
in the feedback gain values and Nissim concluded that a practical flutter 
suppression system would require both leading- and traillng-edge controls. 

The incompressible Pade model of (4,9) is capable of analyzing 
Nissim* s design. The control law of (5.11) is implemented by noting that 
for an oscillating section Ih ^ b/u), giving the control law 


r n 


r n 


"* - 

p 


h 

1 

« 

h 

y 

= l-c,] 


+ — [c,,J 

» 


a 

^ — 

r 

a 


Tlie reference frequency is chosen to be in the vicinity of the flutter 

frequency. For the section of Table V-2, - 76 rad/sec (Fig. V-lc), 

and Fig. V-7 shows the locus of roots as a function of for the 

uncontrolled and the controlled section. (Tho damjnng, Q 0,025 in 

tho figure.) T}io control law stabilizes the bonding and torsion modes 
throughout the range of velocities 0 < U/bw^ < 5.9 but the lending-edge 
ila]! mo(h^ is unstable' tlirf)UghouL this range. The closc'd loop bending mode 
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appi'oitolu's l;lu‘ axis and leads to a sliU. ically diveripait, mode at 

U second case wltli 25 )-ad/st‘C was studietl and a s:lm- 

ilar insi ability was observed in tlio flap mode, lU'iicr, .Vissiiti's desiivn 
pr<Hu>diii'(> is s(>riously deficicniL in nec;lect inn tlie flap de{ci'ees-oj-f la'cdom 

KvaluatiiiK the aerodynamic energy desl<-n t<‘clin i qiK' , the followinf^ 
coiiinu.'nts iwv appropi'ialt' ; 

1. The technique i.s o'/crrly conservative in that it attempts to 
suppress flutter for all ])ossil>le combinations of values of the 
structural and geometrical parameters defining the section. 

2. The technique, which attempts to define a FMC system valid for 
all possible combinations of structural parameters, is incapa-^ 
ble of producing a good design for a single trailing-edge control 
since at least one combinat ion of parameters can be found f oi’ 
which the section is uncontrollable (viz., the section of Table 
V-1). 

3. Tne technique addresses flutter suppression without regard to 
the strucuural properties of the section. To assess the flutter 
boundary of the final design, a standard U-g analysis must be 
performed using the final control law. 

4. The aerodynamic energy eigenvalues would appear to have no di- 
rect relationshi}) to the locus of roots in the s-plane. Thus 
they offer little guidance in design modifications. 

It would ap])ear possible to modif}^ the aerodynamic energy design technique 
to handle the problem of leading-edge flap instability by Including the 
flap modes in the design. The main problem in this extension would bo 
the c<miplexity of optimizing the control laws of a large order system 
ovei’ a large range of k. 

The technique was extended to the design of a flutter sumn'cssion 
system for a delta-wing wind-tunnel model IT)]. Nine flexible modes were 
included in tin.' design but llie leadi!\g- and trailing-edge control surface 
iiKxles were not indudt'd. It Is significant tliat for this largo order 
system it was not possible l.n acdiic'vc' a design liaving all of the eigenvalue 
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oi l positive ovor tho (losirc'd ran^o of k, A1 thour;h t lu' nonpos i 1. i ve 

(1(‘ ti ui tones s of would seen) to imply that fluttc-i’ su|>pi’ess i on hatl 

not been aehlt'vcd, IJ-g analysis verified an increase in the flutter speed. 

The resulting control laws were tested on the wind tunnel model at 
.M O.B, 0,7, O.H, and 0.9, At the first three Alacli numbers, tlie system 
could not be evaluated due to a severe leading-edge control surface in- 
stability. The instability was not encountered at M ~ 0,9 and the 
flutter suppression system demonstrated a significant increase in the 
flutter dynamic pressure. It is suspected that the leading-edge surface 
instability is of a similar nature t- that analyzed in Fig. V-7. 

E. FINITE STATE MODELS OF THE RATIONAL PORTION OF 
AEROELASTIC- SYSTEMS 


The complexity of current design practices and the difficulties 
experienced in implementing designs emphasize the need for simpler tech- 
niques in active aeroelastic control. The ability to calculate unsteady 
aerodynamic loads for arbitrary values of s, coupled with the insight 
gained in the study of the Laplace inversion integral for two-dimensional 
flow point to a new technique of aeroelastic system modeling. This 
technique is acveloped in this section and applied to the study of active 
flutter control of a two-dimensional section. 

The transformed equation describing the aeroelastic system is given 
by (3.1) 

G(s)X(s) :r Gll(s) . (5.13) 

If rigid two-dimensional sections are being considered, G(s) is an nxm 
matrix whose elements may contain nonrational transforms (c.g. , C(s) in 
incompressible flow), while ^Us) is an m-dimens ional input. For flexible 
wings, (5,13) may represent an infinite dimensional operator relation de- 
fined on a Hilbert space. More commonly, the infinite discrete spectrum 
of such a wing is truncated to the. first n modes, 'thus, both cases 
may be treated by the finite dimensional version fo (5,13),.... 'the two- 
d i (IK ‘ ns ional sectiem in L ncomp res s 11)1 e flow will be used to illusti’ate the 
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dc'vol opim>nls of Lliis soc^tion. Tlu‘ .solution for Liu? eompononth of X(^i) 
is 


Xj_(s) 


1 / X 

S N (s)u (s) 
j=:l ^ 


D(s) 


^ - ^ ^ } • • • • n 




where N^(s) is the numerator transik^rm obtained via Cramer rule for 
the i^^^ state due to the input. The inversion integral may be used to 
obtain 


m I n 






) 

CO 

X 

{ 

f - 

-j/ i-n, 

•rt > 

\ Im 

u (re ; 

e dr ( 

J 

o 

.^k 

1 



) 


(5.15) 


The summation inside the brackets has been termed the 'rational' portion 
of the response and is due to the isolated poles introduced by the struc- 
tural degrees-of-freedom while the integral has been termed the ’ nonrational ’ 
portion of the response. 

The examples of the previous chapters have shown that the oscillatory 
transient response typifying flutter phenomena is due entirely to the ra- 
tional portion. The response of the nonrational portion is nonoscillatory 
and decays monotondcally to zero. Moreover, it is a small fraction of the 
total response. In addition, much of the analytical difficulty in study- 
ing the response of the system is caused by this nonrational portion. 
Therefore, a model of the rational portion of the system would be desir- 
able, since it would apparently describe the principal characteristics 
of tlie oscillatory response. 

Tlie method to be described below is similar in spirit to that out- 
lined by Wang [50], who showed that a class of linear distributed systems 
with purely discrete spectra and a finite number of unstable modes could 
ho stabilized with a finite-dimensional linear feedback. However, Wang 
did not address the problem of the construction of a simpler approximate 
model describing the instability. The realization of the model of (he 
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rational portion can bo stated as a theorem. 


TIIKOltKM ; The linear system represented by 


a(s)x(s) - G(s)ll(s) 


(5. HO 


where X(s) is n-dimensxonal , ll(s) is m-dimensional , and S(s) and 
G(s) may contain nonrat ional functions of s, may be approximated 
by the linear, constant-coefficient system 




X = FX + G^u 


(5.17) 


with X e an d u g e"‘. The system matrices are given uniquely by 


,m 


F = TAT"^ 

A = diag(s^, s^, s^) 

Ai ■ 


N 

G V A. 

^ x^l ^ 


where s^, i =1, 2, 
and the elements of the NXm matrices A. are 


N are the isolated singularities of 5(s) 


°jk * '"*jk 


the residue at the pole of the state due to the input). 

Note that the finite-dimensional syste?n of (5,17) represents the 
rational portion of the system of (5.16) which is due to the isolated 
poles at i - 1, 2, N, The following proof assumes that N r. 2n 

where n equals the number of degrees-of-f reedom of the system. It is 
possible that the nonrat ional functions contained in C!(s) may introduce 
singularities into the spectrum of d over and above the singularities 
due to the dimcnision of X(s) (c,f,, the two-dimensional supersonic case 


f ' 
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of Soct, III-A— 2), lu this ease, 2n of the system polos must he seleel.c'd 
to construct the appi’oximatc system. A natural choice would be those 
poles associated with the n structural dcKroos-of-f rcedom. 

The proof of the theoi’em will be constructive, giving an algorithm 
for the construction of F and The result of the following lemma 

will be needed. 

lemma : The matrices of residues, A^, i “ 1, 2, N have rank 1. 

PROOF: This may be seen by noting that (5„15) may bo transformed and 

reorganized as 


N A. 

X(s)-. = ^ • (5.18) 

i=l "i 


with R representing the nonrational portion. The elements of the residue 

matrices, A., are given by N^(s . )/ d * ( s . ) . 

1 K i i 


Now consider the linear system obtained from (5.16) by evaluating the 
nonrational functions contained in fl(s) and G(s) at the pole located 
at s = s^. Denote the resulting matrices as a(s,s^) and G(s,s.). The 
solution of 

fl(s,s^)x(s) = G(s,s^)'U(s) (5.19) 


may be written as 





(5.20) 


where A^ N|^(s^)/D* (s^) are nXm residue matrices. Gilbert [82] proves 
that the A. matrices in (5.20) have rank 1 as is also evident from (5.9) 
At the polo, s =: s^, the systems of (5.16) and (5.19) satisfy the follow- 
ing relations 


-12:^- 


✓ 


Tile ref cu’O 


r>(.S.) 


N;|(«i) . 




nJ(s^) 

D’(s^) 

D'(s^) 

D'(s^) 


showing- that and differ only by a complex multiplier. Therefore 

tile Nxm mati'ices A^ in (5,18) also have rank 1, 

To construct the matrices F and of (5,17) consider the diagonaliz.- 

ing transf ormation X Ty, (The eigenvalues of F are assumed to be 
distinct,) Equation (5,17) is transformed to 


y = Ay -I- gu 


(5.21) 


In teinns of X 


X = TAt'^X + Tgu 

showing that F TAT ^ and -- Tg. The proof will be complete if 

unique matrices T and g can be found yielding these relations. Trans- 
forming (5,21) gives 




(si-/.) ^gu 




S-S 


TIh it’ ( n ‘ 


- 1 24 - 


X(s) Ty(s) -r 



X(s) = 


N 

£ 


1,1 =-=i 


A. = T .g. 

1 •! 1 < 


This shows that it is possible to construct a unique realization (5.17) 
if and only if the matrices of residues can be constructed as the 

outer produce of two vectors. In other words, Rank (A^)= 1. But the 
above Lemma proves that this is the case for the system of (5.16). We 

are free to choose the form of g. as 

1 ♦ 


'l. ' “2 


m 


Then the columns of T are given by (5.22) as the first columns of the 


matrices, A. . That is 

' 1 


T =: 


I t I 

T f T «... I T 
• 1 I .2 » I -N 


where T . - [ A . ] , 
• 1 1*1 


From the discussion following (5,21) 


N N 

■‘■b = £ = £ "i. • 


i--:l 


i-1 


This completes the proof of the theorem. 


It is interesting to note that the amount of information available 

iihout tlic I’ational pc^rtion of (5,1(5) is sufficient to uniquely determine 
2 

the N f Nm unknown elements of F and G Tliere are N values of s 

1 i ’ 


- 125 -- 


ORIGINAL PAGE W 

OF POOR QUALVn 



i ntic^prnilpn i. rlcMiuMils of Uio inatricea oi' rosiduc's A., and N(ui’^i) 


L‘tms l.aul : 


^ f ^ f •••! ,1 


lUuliuu'd ordt.*r approximate models of the systcan of (h^Ui) may l)o con- 
sti’ucted by deleting selected poles since the theorem is usually true 
foi* the case in which the dimension of X is N < 2n. 'this follows 
sinct^ Liie pi- oof of the tlieorem depends only upon the rank of the 
matrices beinp; unity. For the case N < 2n, these matrices will be 
submatrices oI those (‘onsidored in the theorem and will have rank less 
than or equal to one ^ Disregarding the very unlikely occurrence of a 
i‘ank zero submatrix, the above statement follows. 

The modelvS constructed from the algorithm given by the Theorem will 
be termed M*ational r. odels * . The matrices F and describing the 

rational model ot the three DOF section of Table III-3 in incompressible 
flow are given in Table V-4 L'or U'b'..,, 2.9, 


Table V-4 

HATIONAI. MOi)FL FOl^ A TimEF-DEGREE-OF- FREEDOM SECTION 
(M : 0, D. b^o - 2,9, F and in sec ) 

' a ' 1 




The matrices of the Parle model are ^Ivor. in Table V-5 for Uie same case 


Table V-5 

PADE MODEL EOU a THREE DEGREE-OF-FREEDOM SECTION 
(M - 0, “■ F and in sec""^) 



1,0 

0 

0 

0 

1.0 

0 

0 

0 

• 

L 

-9.267 

-10.47 

-0.9598 

12.30 

-15.52 

-4,820 

-25.61 

16.50 

-8.755 

0 

0 

0 

1.0 

0.9 

0.1487 


0 

0 

0 

-583.2 

774.2 

-1611 

1.0 

- 100.2 


0 

0 

_0 

355.3 

•8939 

115712 

0 

0 - 


Table V-6 compares the transfer functions derived from the rational 
model with those of the Pade model. Since the F matrix of the rational 
model is constructed by performing a similarity trunsfoi'mation (5,21) on 
the matrix of the exact eigenvalues of (5.16), the rational model repro- 
duces these exact poles while the Pade model gives a good ap|)roxiinat ion 
to these poles. In addition, the Pade model also contains the two aug- 
mented state poles^. The nature.- in--which the residues used to construct 
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Table V-6 

COMPARISON OF TRANSFKU FUNCTIONS OF RATIONAL AND l■Ab^'; MODELS 
(M r.: 0. u/bei r: 2.9, polos and '/.oroos in see ^ ) 


Rational Model 


Pndo Model 


bending 

-3.711 + 

168,252 

-2. 

,659 

+ 

169, 

149 

torsion 

-15.049 + 

180.171 

-17. 

,485 

± 

179, 

122 

flap 

-3,431 + 

1340,38 

-3, 

,414 

+ 

1340 

.26 




-12, 

,523 







-74, 

,099 





2) ZEROES 


a) Rational 


Model 

Gain 


7.2785 


-95.518 


-3.3904 

x^/P„ 

-8648.6 


(-4.09±i81.59) (-65.3+1364.68) (+107,43) 
(-3.47+180.18) (+110.91+11173.1) (-30. 70) 
(+2.82+171.86) (-72.0iil51.21) (-2457) 
(+3.3 +171.55) (-56.48+1152, 21) (+10,44) 


b) Pade 


" . - . ... — ■ ■ 1 „ - ■ .... j 

Mode 1 

Gain 



355.26 

(-3.25+180. 28) (+501 . 72) (-334.66) (-220. 3) (-15.53) ' 

a/p 

-8938.89 

(+3.84±i72.37)(-59.84±U62.29)(-56.09)(-12.01) 


V and were evaluated places constraints on the system realization ^iven 

by the rational model which are seen in the structure of the submatrices 
of F iti Table V-4 , Th<‘ equation for Is 


X. 4 7,279 


(5.23) 


/ 







I 


liullcnting thc> cxpectcMl rcilatloa botwocn the utatos mode! inn l> ( •'j ) and 
li The small ' Icedforward ’ tenn, 7,a79f’'^, c;au;:es the ntvlat ion 

Xj f. and is tlue to the unmodolod nom’ational portion oI (5.]r>), This 

effect is also evident in the zeroes of Table V-6, The Padc^ model zerot's 
satisfy the relation (h/f\)(s). i- s(h/p )(s) while (X /['• )(s) / 

«(Xi/P^)(s), The magnitude of the feedforward terms of the upper sub- 
matrix of of Table V-4 is direct )y proportional to the relative 

magnitude of the nonrational portion to the rational portion of the re- 
sponse. It is interesting to note that the real zeroes of the rational 
model indicate phase changes of approximately 90” between and , 
and X and X respectively. 

The usefulness of the rational model must be evaluated by its ability 
to predict, the response of the system in the bandwidth of interest; that 
is, at frequencies near the flutter frequency. Table V-4 shows that the 
zeroes near the axis in the vicinity of the flutter frequency (u) = 

70 rad/sec) agree well between the rational model and the Fade model. 

Frequency responses of the rational model are compared to the exact model 
in Fig, V-8. The rational model agrees very well with the exact model 
in the frequency range of flutter, with the agreement deteriorating \vith 
increasing distance from the system poles. From this comparison and that 
of Fig, IV-4, it is concluded that the rational model and the Fade model 
are both capable of predicting system response at frequencies near the 
flutter frequency. 

In closing this section, it should be noted that the rational model 
is not restricted to the two-dimensional incompress ible flow case. It is 
equally valid for compressible three-dimensional flow when used with 
truncated normal mode structural representations and aerodynamic loads 
calc^ilated for arbitrary s. The advantage of the rational model that 
il does not require augiaented states whereas the Fade model does. Also, 
the rational model might be expected to give better performance than the 
Ihido model for points well removed from the b. axis. 

f 
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Exact model 
Rational model 



(a) PLUNGK 

FIG. V-8 


COMPAllKSON OF FRFQUKNCY RKSPONSES OF PLUNGK AND 
1'ORSION DUlv TO FLAP DEFLECTION OirJ'AINEl^ FROM THE 
I'XACT MODEL AND FROM THE RATIONAI. MODEL 
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CCjNTHOJ, CJK A}yUOKJ_AaTIC,..iiYSTJ’:MS 

la this st.»ct;i.ou actlvt' i'lutUar (losi.|^’n lochniqiu's will l>o stiulled 
usLni;' Ihv ‘•alional modi'l ami the* l^ulr modo I . Although i:l )0 spec 11 ic prohlcan 
undiM* eonsidt>raLion will )>o ihc coiitzv)! o£ a I’lultcr mode*, the sauio looh- 
niqiu's aiH> ajjplLcfibio to any of thr couLrol confl^urod voiiiclo (CCV) con-^ 
c('pLs addressiun* dynamic Htructural or aircraft response. 

Several of the most promijicnt examples of noroelastic control have 
been designed using augmented state Pade models [2], [3], which resulted 
in (juite large order systems. Attempts to apply modern optimal control 
techniques to these models have not had great success due to the require- 
ment of feeding back all of tiie states of the model. The matrix Pade 
ap})roximant method of Vepa promises to alleviate tliis problem somewhat 
by greatly reducing the number of augmented states. The rational model 
holds further promise in that augmented states are completely eliminated. 
The cost of this advantage is a certain ambiguity in the relationship 
of the states of the model to the physical measurements of the structure, 
ITiis ]>roblem may be assessec^. by checking the performance predicteci by the 
rational model with the corresponding performance obtained with the exact 
model using the same control law. 

In -the two-dimensional incompressible case, it is possible to compare 
the two models since exact unsteady airloads are available for arbitrary 
s (vSec, II-D) and the Pade model of Sect, IV-A involves only two aug- 
mented states. In both oases, the model is given as 

X = FX G^u • (5,24) 

An appro))riato performance index foi the flutter problem is that of the 
optiiiial )'(‘gulator [37] 


CO 

1 r T’ T 

4 \ I X AX I u Hu dl 

^ •’o 


^ 5 , 25 ) 


M i n imi /.a t ion ol J satisfying Ihe c:onslt'aini of (.T ^ 24 ) is achieved via 
tin? s lead y~s L a !(' solution of tlio matri;. Kiccati ecpiatlon 
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,S -,SK - F^s I SG^H"*(i|S - A 

wlu'i’o S(0) 0. 'Plu' linear roodbacM; control law Is n‘tv(‘n in terms r)t 

o 

the steady-state Ui(*eatl matrix, S , as 

u(l) : -b''‘'g'J’s"x( 0 - cx(i') . (5.27) 

The control weighting matrix B must be positive definite, while the 
state wei{.^hting matrix A must bo positive scmidcf ini to , Hall and Bryson 
[84] describe a digital compute!* program well suited to the solution of 
this problem. This program, OPTSYS , was utilized in the control law de- 
signs of this section. 

The i^roblem of choosing the weighting matrices A and B remains, A 
basic result of the theory is that if A rr o (i.e,, zero-state weighting), 
the action of the resulting control law upon the closed loop eigenvalues 
is to leave unchanged all stable eigenvalues, while open loop unstable 
eigenvalues located at s = u + iu) are reflected about the ito axis to 
s r: -UL 4- iu. For structures with slightly supercritical flutter modes, 
this zero-state weighting technique is an attractive design method since 
the result is a modestly stable controlled mode. However, for larger 
values of the supercritical flutter speed, the technique leads to unrea- 
sonably high damping of the flutter mode. Also, the method cannot be used 
to improve the damping of subcritical flutter modes. 

Anderson and Moore [85] describe a method which can be used in con- 
nection with zero-state weighting to restrict the maximum value of the 
real part of all eigenvalues. Equation (5.24) is transformed and the 
change of variabl(\s s* n s -f v made giving 

Ls'I - (FI-vI) ]x(s' ) == G^'ll(s) . (5.28) 

The result of the t rans format ion is to shift all eigenvalues of F v units 
in the direction of tlic positive real axis or, equivalently, to sJiift the 
origin v units in the direction of the negative real axis. K(iuation 
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i lul i t'.a t.liat this (M ^oi^valuo .sJU.l’t can bo aocoiii])] i sIumI by a^kl i ng 
V Lo Lilt! (li;it;oii!il cli'iiitnls of K. llonce, if F contniuH an unsLal)l(> 
o Ini'iiva liii' at s |i ^ iw and II, in dt'Mi red to constrain al l cin<*nval uc.s , 
such that Ui’(X) < 1], this will he accom|ilisli<'(l |jy the o))Limai reg- 

ulator .solution witli /ei'o-state weinhtinf’' foi/ the system 


X :: (F ( vl)XG^U 


(h. '/!!») 


with V -e(lt-Tl). Fin'envalues , witli ReO,^) < T], arc unaffected 

by the resultlni^ control law. 

The optimal regulator solution was obtained for the section of Table 
III-3 at - 3,23 using the rational model and the Fade model, (The 

flap chord, c - 0,5, for this example.) Figure IV— 3 compares the open 
loop roots of the two models., and shows that the section damping is 
Q ~ -0,09. The weighting matrices were A = 0, 13-1. Table V-7 gives 

the feedback gains and the open and closed loop eigenvalue locations. 


Table V-7 

OFfliVlAL REGULATOR GAINS AND EIGENVALUES FOR A THREE DOF SECTION 

(M = 0, U/b(jjQ, = 3, 25, A = 0, B = 1, poles in rad/sec ) 

Case I: Rational Model 


C= 12.901 -2.197 -0,09393 i 0.04124 0.007558 0.00010641 




bending (flutter) 


tors ton 


Case II; Fade Model 


Open Loop 


+ 6.420 171.03 
-28.68 ± 173.56 
-14.59 ± 1339.9 


Exact Closed Loop 


-4.975+ 169,94 
-28.68 + 173.55 
-14.59 ± 1339.9 


C :: [2.517 -2.519 -0.1076| 


l0.0450 0,01168 0.00051831-5.778 0.1436] 


boHdinu ( J lutter) 

tors ion 

flap 


Open Loop \ 


+6.987 » i71,01 
-31.27 » 172.82 

-14.53 1339.6 

-1 


PJxact Closed Loop \. 


171.48 

i69,45 

1339.7 



ORIGINAL 
OF POOR QUAIvITY 










i 


Tiu> rods ‘ opcsi loop» in (lu‘ Tsl)Io arf> (h<' rsi^i -nva 1 ui‘s of thi* 

appropriate unenntrollfsl K matrix (vi/.., Kq. 5,1V TfH’ tiio I’alional 
modol, and Kt{, 4,9 ior thci Pad^^ model), Tlic roots tormc'd ’exact clf)sc‘d 
loop’ wore f)btai nod by i mploinont Inp,' tin* feedback control laws in the exact 
system equations and locating* tlie exact closed loop roots by iteration 
as described in Sec, For the rational model, (5.15) becomes 


|,G(s) -GCiX(s) = Gll^(s) (5.30) 

where ll(.s) ll^(s) + CX(s). To implement the Padc< model control law, 
estimates of the two augmented states are required. This was accomplished 
by adding to (5.13) the known structure of the augmented states. In 
tile notation of (4.8) this yeilds 


’G(s)-G(C^+C2s) { -GC, 


X(s)‘ 


Ei+EgS 


P J 


|Xp(s)_ 



(5.31) 


where C , and C contain the gains associated with x, x, and x . 

1 . j - - -p 

If the linear models represented by X = FX + G^^u exactly described 
the dynamics of the section, then the closed loop roots would be Identical 
with the open loop roots except for a sign change in the real part of the 
unstc^ble roots. Deviations of the roots away from this condition indicate 
tlie presence of unmodelcd (nonrational ) effects. 


Table V-7 shows that feedback gains from corresponding states of the 
two models are comparable and both designs stabilize the flutter mode. 
Except for the h feedback gain, the magnitudes of tlic gains of the 
rational model are less than those of the Pad^ model. The distance from 
the exact closed loop pole location to the predicted pole location is 
an indication of the ability of the models to deal with the unniodeled 
portion of the system. These distances for the flutter mode are; rational 
model, /vS - 1.56 rad/sec; Pade model. As = 1,29 rad/sec. The corres- 
ponding distances for the other two modes illustrate a basic difference 
between tlie models. The ratii.nal model is an exact model of tlie rational 
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portion oi' the .system at a particular condition and the regulator 

design with /ero-.statc weighting leaves the exact stable ])oIes unchanged. 

The I’adc? model is an attem])t to approximate the system throughout a given 
bandwidth (region of the s-plane), 'fhus the closed loop torsion and 11a]) 
modes are displaced 3.74 rad/sec and 0.71 rad/sec respectively from their 
predicted locations. This effect of Fade models may be of concern in tlie 
design of flutter suppression systems for multi-mode structures in which 
there may be several marginally stable modes in addition to a flutter mode , 

The cliuractoristic of the rational model of matching the open loop 
rational portion of the system exactly emphasises the perturbation nature 
of control laws based upon this model. This implies that the deviation 
between the predicted and actual root locations will increase as the 
distance by which the flutter mode is moved increases. The deviation 
indicated in Table V— 7 would seem to be acceptable. If the deviation were 
unacceptable, a second rational model could be constructed for the S 5 \stem 
resulting from the use of the first control law and a second regulator 
design performed, giving a second control law. If this attempt resulted 
in a satisfactory design, the final control law would be formed by the 
sum of the two control laws. Thus the rational model can be used in an 
iterative fashion, whereas a corresponding capability is not apparent in 
Fade models. 

Figure V-9 indicates the effect on the open loop eigenvalues of the 
Fade model of incorporating the Fade gains of Table V-6 one at a time. It 
shows that the main cont ributors to the s tabilizat ion of the flutter mode 

are the h, a, and h gains. The h and p gains destabilize the flap mode, 

0 

an effect which is counterbalanced by the Ct, P, and Ot gains. Inter- 
estingly the augmented state gains have little effect on the flutter mode, 
but they do influence the remaining modes. 

Figure V-10 shows the effect on the exact closed loop pole locations 
for off-deslgn airs])eeds from U/boj^ 0.5 to 3.75, Both closed loop sys- 
tems are Unstable below the open loop flutter speed of 

fleeting the near uncontrollability of tlil.s section by tlie Irai ling-edge ^ 
control surface at this airspeed. Above this airspeed, both control 
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laws pi'ovldc flutter mode control until a static divergence occurs at 
approximately - 3.80. The Fade model maintains slightly better 

stability than the rational model throughout this airspeed region. Tliis 
is due to the ability of the augmented model to approximate the nonrat tonal 
effects over a wider bandwidth than the unaugmented model. However, this 
capability requires the complexity of augmented states, with the attendant 
problems of state estimation. It should also be noted that although the 
rational model control law was designed utilizing a 'nonphysical' model 
(viz., Eq. 5.17), the performance indicated by the 'exact closed loop' 
pole locations was obtained using measurements of real physical states 
(h, a, p, h, Q:, and indicates the performance of the system under 

actual operating conditions , 

The optxmal control of the three DOF section was also investigated 
in supersonic flow. The section of Table III-4 was studied with the 
rational model and the matrix Fade model of (4.19), Table V-8 gives the 
optimal gains and eigenvalues for the two designs. The weighting matrices 
were A - 0, B = 1. Again, the corresponding gains of the two models are 
comparable with the rational model gains having smaller magnitude in all 
but two cases. For this compressible Pad^ model, three augmented states 
are required; one for each degree-of-f reedom. The exact closed loop 
poles were not calculated for the Fade model. The agreement between the 
open and closed loop poles of the rational model indicates that the un- 
modeled effects are slight at this Mach number. Figure V-11 shows the 
effect of off-nominal values of M on the closed loop poles when the 
M " 2 feedback gains of the rational model are held constant. The 
figure indicates that the flutter Mach number has been increased from 
M ^ 1.8 to M = 2.2. Comparison with Fig. IV-6 indicates that the control 
law also stabilizes the flap mode at the lower Mach numbers. 

As a final design case, the four DOF section of Table V-2 was 
analyzed using the rational model. Figure V-lc illustrates the nature 
of the flutter mode which was studied at U/ bw 3.55. The flutter 
mode is unstable with a damping of Q i -0.23 and the airspeed is 25 per- 
cent above the flutter s))oed. Also, this section is nearly uncontrollable 

by the trai 1 ing-edge control surface at the flutter speed of U/bu) 2.84. 

' a 
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Tabic V-8 

OPriMAT, RKGULA'l'OR GAINS AND KIGENVALUliS FOR A THID'JE-DEGREE-OF 
FREEDOM SECTION IN SUPERSONIC FLOW 
(M r- 2,0, A = 0, B r 1, poles in rad/sec) 


Case I: Rational Model 

C = [0.300 0.132 0.00401 j -0.00513 0 

1.00932 . 0.000349] 

Mode 

Open Loop 

Exact Closed Loop X. 

flutter 

-4.036 ± 174.67 

-4.078 ± 174.24 

bend ing-tors ion 

-15.69 ± 170.67 

-15.69 ± 170.67 

flap 

-4.353 ± 1372.6 

-4.346 ± 1372.6 


Case II: Fade Model 

C = [0.46 -0.046 — 0.078 j 

I -0.004 0.010 0.0004 ! 0.034 0.086 -0.032] 


Mode 

Open Loop 

Exact Closed Loop 

flutter 

bending- tors ion 
f lap 

-4.249 + i74.67 
-15.86 ± 170.60 

-4.374 ± 1372.8 
-482.8 ± 13323 

-1855 

— 


lable V-8 gives three designs accomplished wltli zero-state weighting. 

For Case 1. B diag (1,1), weighting the leading- and trail j ng-edge 
control motions equally. Tlu; deviation of the closed loop llutter mode 
from its anticipated location (.<- -16.66 168.08 i-ad,sec) is accept- 

able but the gains relating to the leading-edge control arc significantly 
higher than tliose associated with the trai ling-edge control. The regulator 
solution lias designed a control law calling for inort' mot ion by the leading- 
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THE EFFECT OF OFF-NOMINAL VALUES OF M ON THE CLOSED LOOP POLES OF A 
THREE DEGREE-OF-FREEDOM SECTION IN SUPERSONIC FLOW USING THE 
OPTIMAL REGULATOR GAINS OF TABLE 5.81 


I 


c‘.oi\U’()i .sinci' tliiK sui’fiuu' producos urcMitor loads X'or a nivon dcXloc- 
t. Ion tlum Llw' tra 1 1 .lnn-(Hl(j;c control. Tills is undos i rabli^ since' tlio 
lcadin^;-od^•c control is also sut),)oct to muedi n'roator liin(;(' momonts and 
would consuino o eir res pond inu'ly greater power. 

Tims the weighting on the leading-edge control motion was .increaseel 
and Case II was designed with H ^ dlag (1, 1(5), This change implies 
that tliat leading— edge control motion is l(i times more 'expensive' than 
tra iling— edge conti’ol motion. Table V— 9 shows that tlic deviation of the 
closed loop flutter pole from its expected value has almost doubled its 
value in Case I but is still acceptable. Also, the trailing-edge flap 
gains are roughly doubled from those of Case I while the leading— edge 
flap gains have been reduced by a factor of 5 to 10. The increased 
weighting on the leading-edge control has produced a design requiring 
more activity by the trailing-edge control. The increased deviation of 
the flutter mode from its expected location is probably due to the prox- 
imity of zeroes associated with the trailing-edge control. 

The final case of Table V-9 illustrates the Anders on-Moore technique 
of axis shift in the s-plane to achieve a specified degree of stability. 

The desired damping sought by the regulator solution for Cases I and II 
is quite high (Q = 0.23) and a significantly smaller damping would be 
acceptable for high frequency structural modes. Using the Anderson-Moore 
technique with v —5 rad/sec will result in a control law which attempts 
to place the flutter mode at s r; —0,06 + 168,08 with a damping of Q - 
0.10, Case III gives the resulting design for B - diag (1, 1), For 
thi.s more modest design goal, tlie deviation of the exact closed loop pole 
from its anticipated location is only 1 rad/sec, les than half of the 
tleviation of Case I. 

Figure V-12 shows the migration of the bending and torsion modes 
for Cases 11 and III as a function of U/boi^. The feedbexek gains of 
Table V-9 were held constant and the exact closed loop eigenvalues located 
bv iteration. The design of Case II is unstable at low values of U/bi.:^ 
and gain scheduling would bi' reciuiri'd to achit'Ve acceptable periormance, 
Kviih.'Htly, |•l•ound checkout of this fluttci- suiipiu'ssi on system would bi' 
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rablo V-9 

OrriMAL ItKGUI^TOU (JAIN.S AND KIOKNVALUliS FOR A FOUR-DEGRRl 
FUEEDOM section using TICE RATl'ONAI. MODEI- 
(M f). A 0. U/b,„^ . 3^56, poles m rnci/soc) 


Modt 

1, bond inf; (flutter) 

2, torsion 

3, L.£« fUp 

4, T,E, flnp 


Open l<oop 
416, G6 1 IG8.08 
•*33,62 ± 163, C3 
-i4,7l i 1254,38 
-02.76 ♦ 1605.10 


Case I: 0 a dlag (l.l) 


0.62 0,173 0.0105 -0.0760 j 

0.787 0,316 0,0360 -0,139 i 

j 0,00333 0.00823 0.000611 -0.000987 

I 0.000586 0,0115 0,000730 -0.00148 


Wode 

U bendlnc (flutter) 

2, torsion 

3. L,E. flap 
1. T.E. flap 


Exact Closed Loop 
-15,29 i 169,77 
-33.52 t 163,77 
-14.71 ± 1254,38 
-62,75 t 1605.10 


' ' ' 

C«s« 11: 

B = dlagd, 16.) 





ri.479 U.451 

0.0282 

-0,187 1 


C 

* 


I 



Lo.118 0.0746 

0.00531 

-0,0204 I 



j 0,0075 

0.0198 

0.00124 

-0.00239 


1 0.000121 0.00171 

0.000109 

-0.00022 

Mode 


Exact Closed Loou 


1, bending 

(flutter) 

-12,92 t 

166,08 


2, torsion 


-33,62 t 

163.63 


3. L.E. flap. 

-14,71 t 

1234,38 


4. T.E, flap 

-62.75 ± 

1605.10 



Case III: b * dlag(l. 1.), 

V « -5 rad/scc 



ro.439 

0,0899 

0.00485 

-0.0494 



C e 





10.649 

0,3U 

0. 0219 

-0.0899 



1 

1 

0,00277 

0,00574 

0.000354 

-0.0006771 


1 

0,00101 

0,00785 

0,00501 

-0. 000997 J 


Mode 


E.xc.>st ^Closed Loop ) . 


1, 

bending (flutter) 

-6,079 

t 168.91 


2, 

torsion 


-33.62 

t 163.63 


3. 

L.E, flap 


-U.71 

t 1254,38 


4, 

T.t. fUp 


-62,76 

± 1605.10 




0 


(a) A - 0 . B = DIAG (1 , 16) 

FIG . V - 12 THE EFFECT OF OFF-NOMINAL VALUES OF U/bco^ ON THE 

CLOSED LOOP POLES OF A FOUR DEGREE -OF -FREEDOM 
SECTION IN INCOMPRESSIBLE FLOW USING THE OPTIMAL 
REGULATOR GAINS OF TABLE 5 .911 , HI 




Torsion mode-\ 



FIG. V 12 


CONCLIIDFD 



(li.rric'ull (luH to llK' Liuwlia cnupiiiif;' nl* tin? modt*.s . V~l.^a also 

shows Lho inriiuMioo ol‘ Lh<* unooutroliablo modo i ti llu' sonsiLivtiy ol tlu* 
locus near the o))ou locjp riut:tcr volr>city (U/b(.;^ 2,84), The porrori'j- 

ance oi’ the ccaitrol law oi' (biS(^ III at of r-iiomiiial values of U t)(.^ 

a 

(Klk, V-lil)) Lndlcuites a smooth variation of the' locus tlu’f)UKhout the 
v<*iooity ranpa* studied, ihis is due to the inci’eased autlioi’lty allowed 
tlu' lead iuj;’- edge fla]> by the reduction of the weighting on y for Case 
III, 

Both the torsion and the plunge modes of Cases II and III are stable 
for airspeeds well above U/b(.) : 4,0, However, at U/bw - 3,98 the 
section becomes statically divergent due to the emergence of a real 
]jositiv(* root. As the airspeed increases beyond this divergence speedy 
this root becomes more positive. This static divergence of the actively 
controlled section is of the same nature as the static divergence of tlie 
uncontrolled section studied in Sect. IV-D, Wliereas divergence of un- 
controlled sections usually occurs at iiigher airspeeds than flutter, it 
is seen that active control of flutter may reverse this condition. Hence 
tile behavior of active flutter control techniques should always be in- 
vestigated at the zero frequency condition where static divergence ocv'urs 

The examples given in this chai^ter illustrate the application of 
modern optimal control theory to the design of flutter suppression system 
The matrix Bade model and the rational model are botli capable of jiredict- 
ing closed looj) perf oia^nince , The disadvantage of the l^ade model is in 
tlie estimation of the a\igmented states which would be required to imple- 
ment the control law. The use of both loading- and t railing-edge control 
surfaces will obviously simplify the problem of stabilizing the flutter 
mode but the additional control surface introduces other problems of 
stability and i30wer requirements. It should be noted that the nominal 
i ticomprc'ss ible flow S('ctiot\ I lives ti gated in this chapter re? presents a 
W(3i’s t c ase dc?s ign situation, in that it was near 1 y uncont reliable by the 
t r*ai 1 i ng-edge control suj’facu/ at tlie fluttei* velocity. In a realistic 
design situation in wliic.h a flutter suppression system is to bt' desigmnl 
tor a paiMMcular tliglil, condit. i(3n f>r range conditions, i 1. may well b(> 
possible to locate ainl size a single t ra 1 1 i ng-(‘(lgo c'ontrol surihun* to 
achieve t.lu' dfsign goal. 
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Chupl.cM- VI 

SUMMARY AND RKCOIvIMENDATIONS 


ORIGWAl. PAGW ^ 
. )K POO?s QUATiW 


A* ^U^l^lARY 

1 , The transfer funct Lon relating airfoil motions to the airloads 
due to circulation in t\vo~d irnensional incompressil^lo flow Is derived and 
is identified as the generalized Tlieodorsen function , valid for arbitrary 
airfoil motions , 

2, Examples of exact airloads due to transient, stable airfoil 
motions in two-dimensional incompressible flow are p:iven, 

3, It is shown that the solution of the unsteady aerodynamic par- 
tial differential equation for compressible flow contains a portion 
which is linear with respect to the transformed ai-rfoil motions and a 
portion which is linear with respect to the initial conditions of these 
motions. The stability or flutter problem is solely dependent upon the 
first portion which is described by a partial differential equation 
formally identical to that of simple harmonic motion with the replacement 
of i(A. by s. It is conjectured that computer programs which calculate 
simple harmonic airloads may be modified* in a fairly straightforward 
manner to yield generalized airloads. The conjecture is shown to be 
true in two-dimensional supersonic flow and the derivation of generalized 
airloads f Ox ^ his case -is— given, 

4 , 1 h<i generalized airloads are incorporated into the equations 
of motion and the exact locus of roots calculated, giving quantitative 
results regarding subcritical and supercritical flutter conditions. 

5, Examples of exact airfoil responses due to command inputs are 
given. The responses are shown to be composed of portions due to rational 
and nonrational transforms. It is shown that the oscillatory motions 
typifying flutter phenomena are due entirely to the rational portion of 
the response , 
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'I'lio ability to cal aula to j^LMicnili/.ocl aorod ynamic Itauls allows 
Llic (walualloii o£ appi*oxiinatc tccnni ((U<\s of cal (uila I i up; Idicsf? If)ads, the 
\i ,'i\ Jou('s’ approximation for 1 neompress l.blc flow and the matrix Ibadf^ 
approximants of .supersf)nic I'low arc shown l;o pivc accurate* airloads lh)r 
arbitrary motions well romoveci from tlie Wu axis. 

7, It is sliown that static diver go nee of typical sections in i n- 
compressible flow occurs due to tlic emergence of a real positive pole 
of the system transfer function, 'Phis pole occurs in addition to the 
original structural poles and is also predicted by PadS approximant 
methods if the low frequency behavior of the approximants is valid, 

8, It is shown that the aerodynamic energy design technique for 
flutter suppression, which attempts to define flutter mode control laws 
valid for all possible combinations of structural parameters, has diffi- 
culty treating the typical section with a single trailing-edge control 
surface due to the possibility of the section being uncontrollable for 
some selection of parameters. Also, the technique must be extended to 
include control surface dynamics in order to circumvent a problem of 
leading-edge control surface instability, 

9, A theorem is given stating the possibility of constructing a 
unique finite dimensional, linear model of the rational portion of the 
system response which does not rcqui»"e augmented states. The proof is 
constructive, giving an algorithm for the derivation of the ^rational 
model ’ , 


10. Optimal regulator flutter mode control systems are designed 
using the rational model and the Pade model for incompressible and super- 
sonic flow. Although the rational model represents only a portion of the 
total respotise, it is shown that perturbation feedback control based 
upon this model yields acceptable flutter mode control systems. It is 
also shown that active flutter control techniques may result in systems 
with divcrgeuce speeds Ijelow the actively controlled flutier speed. The 
l;eliavior of sueh systems should always be investigated at the zero fre- 
qiu*ney condition whore static divergence occurs. 
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B, RKCQMMKNDATIONS 



1 , W Liul 1 unm'l s ludios rmd f 1 ight tests should bo per for mod to 
establish the validity of the transient responses presented herein, and 
to Inves t i a‘a it-' ilitJ effects of the I’ational and nonrational portions of 
the response. 

2, Existing’ computer programs which calculate sim])le harmonic 
airloads could be modified to calculate generalized airloads and the 
results compared with existing solutions, experimental wind tunnel, and 
flight data. 

3, Tlie possibility of obtaining approximating functions of gener- 
alized aerodynamic loads over a region of the s-plane should be studied. 
These approximations may be of the form 


m 

q (s,M) = ^ f (0,M)g (u;,M) 

and would be useful in calculating the locus of roots of the system. The 
merits of this exact root locus technique versus traditional U-g flutter 
analysis should be studied. 

4. The relative merits of rational models and Pad^ models for the 
analysis and design of aeroelastic sj^stems require continuing investigat ion. 

5, The Laplace transform techniques used herein may be applied 

t'o the gust problem, leading to a unified theory of the control of aero- 
elastic systems excited by turbulence. The ability of the finite state, 
linear, * rational model* of such systems to predict the main features of 
the total response may serve as a base for future applicati .ms such as 
gust alleviation and vehicle ride control. 

b, 'Plio problem of estimation of tlie states of the r*ati ona 1- model 
from measurements of the physical airfoil requires careful study since 
the states of the rational model do not correspond directly to physical 
nieasurements , To obtain complete cor^^espondence , the nonrational portion 
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(or :m osi lmato llioreol) must be included. Hence, tlic operation oi flllerH 
02’ obh’oi’vers for state j'eeons true t itni based upon the i-ational model mus L 
be carefully evaluated. 


7, riie concept of rational and nonrationnl portions of the aij'foil 
response may lead to improved estimates of flutter mode dampiiiK from 
flight tests, Tf the nonrational portion of the response can bo estimated, 
subtraction of this estimate from the total response measurements would 
provide estimates of the rational portion. Application of parameter 
identification techniques to this portion may give improved damping 


estimates . 



Appendix A 


EQUATIONS OF MOTION 


Thr equations oi motion of the section shown in Fi[^. II-l are de- 
rived from La^rango^ s equations 




where the kinetic energy is 


,3 A 


2 •'-i P 


The airfoil deflection for the Sect, of Fig, II-l is 


-h-(x-a)a -(d-x)yil( d-x) - (x-c)pll(x-c) , 


The potential energy is stored in springs attached et the control 

surface hinge lines ^^id et the elastic axis * 


r ^2*2 *2 *2 *2l 

mb h + i^a + i^p ...+. J + s^bha 


+ Spbhp + S^bhT + [ip + Spb(c-a)]a(3 
+ [S^b(d-a) - I^]6^ 


and the equations of motion tor the section oT Fig. II-l arc 
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mbii I S w, f f- f P (A/^) 

S bVi + I rl; , [I IS b(c-a) |(i fl.S lXd-:i}- 1 ]7+k /i' = M' (A.o) 

u r/ ( > I j r f ^ / 

S bii +[I I S b(i;-i\) jlj; I I I 1<J5 = (A.Y) 

i-. l-s (i (. (. 

S-^bh i-[S^b(d-a) - I^]a 1 I^T i- k^Y = . (A. 8) 


In mechanizations of such typical sections in a wind tunnel or on 
a wing, the control surfaces are commonly controlled by electrohydraulic 
servos as described by Edwards [86] and Bergmann [8?]. Then = ky = 0 
and additional terms giving the hydraulic pressure control torques would 
be added to (A. 7 ) ^^nd (A. 8 ). Edwards [86] derives the equations of such 
a hydraulic control system and gives the transfer function from control 
surface position command to control surface position as 


i-(.) 

c 


/=2 2C.S 

{t. S+ly ~ + ii- + 3. 

P V 2 w. 
h 


(A.9) 


The hydraulic mode, described by ul and C , is. typically a lightly 

n n 

damped mode well above the bandwidth of the servo (given by rad/'sec). 

To retain the control surface dynamics in the equations of motion 
without requiring attention to the servo loop dynamics, the artifice of 
control surface springs will be retained and viscous damping terms will 
be added to the control surface equations to provide stability. Also, 
to provide a mechanism for control surface positioning, the control 
surface spring constants will be multiplied by the difference between 
surface position and commanded surface position. Thus the terms 
and k 7 in (A, 7) and (A, 8) are replaced by k ) i- 

Y pc} j 

k^(y-“7^ ) ^ T G sp oc 1 1 V c 1 y , The selection of • 


^ Vk^/l!^, C o ibe flap dynamics to approximate the 

7 7 7 i‘> ) 

hydraulif: position servo loop dynamics of (A. 9 ) 






The aerodynamic loads acting on the section of Fig, II-l may be 
derived from those given by Theodorsen [ll] and Theodorsen and Garrick 
[51] for the section of Fig. A-1, This section has trail Ing-edge aileron 
and tab control surfaces which are aerodynamically unbalanced. Using 
a superscript bar notation to identify quantities related to the section 
of Fig. A-1^ the coordinates of the two sections are related as 
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1 

0 

0 

b(d-a) 


h 

a 


0 

1 

0 

-1 


a 
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0 

0 

1 
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0 

0 

0 
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Vx 


(A, 10) 


while the loads acting on the sections are related as 


L = 


M 






M 
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b(d-a) 
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-1 


0 0 
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0 1 
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(A. 11) 
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FIG. A-1 


DIAGRAM OF A TYPICAL SECTION WITH 
AERODYNAMICALLY UNBALANCED 
AILERON AND TAB 
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Appendix 13 

UNaTKADY AKltODYNAMIC LOADS FOR 33VO- DIMENSIONAL 
INCOMPRESSIBLE j-LOW 

This appendix summarizes Theodorsen’ s [ll] derivation of the un- 
steady airloads as presented in Bisplinghoff [ 7 ] . 

The airfoil lying between x* = -b and :c* = b, as shown in Pig. 
B- 1 , is mapped onto the circle of radius b /2 by the Joukowski trans- 
formation, 

2 

x-x- + iz^- = (X + iZ) + 2 . 

4(X + iZ) 

The correspondence between points on the airfoil and points on the 
circle is x> = b cos =, z* = 0. Solution of (2.19) subject to the 

boundary condition, (2,13) (which is Laplace's equation in the plane) is 
achieved by superposition of elementary solutions of Laplace’s equation. 
To satisfy the boundary condition, a distribution of sources is placed 
on the upper semicircle and a corresponding distribution of sinks is 
placed along the lower semicircle. The source strength distribution 
required is 


H'^(x*, t) = 2 w^< (x*, t) . 

cl 


(B.2) 


This noncirculatory source-sink distribution gives the tangential vel- 
ocity at the circle as 


nc 


„ w* sin <t)d<l) 

( 0 ,t) = I ( . 

> .1 COS (t>-COS 0 


(B.3) 


The noncirculatory velocity potential on the upper semicircle and the 
pressure dirrorcncc on the nirfoil nro 




FIG. B-1 


CONFORMAL TRANSFORMATION OF THE x*-z* PLANE TO 
THE X-Z PLANE 
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The Inni'ontial velocity, (B.3), cvnluatcd ;U (he t railing- edf;c 
( ' ^ O) is nonzero for general airl'oil motions and Kutta's condition of 

smooth flow off of the trailing-edge is violated. 

To satisfy the Kutta condition, Theodorsen employed a bound vortex 
distribution over the airfoil chord, and a vortex distribution over the 
airfoil wake. Figure D-2 indicates the vortex flow for an isolated 
vortex pair in the X-7, plane. To maintain the circle as a streamline, 
a vortex of strength at X = b /4x is paired with a vortex of 

strong(.!> at X = x* 

Von Knrman and Sears [60] show that the corresponding situation 

in the x* - z^ plane consists of a vortex of strength -F at 

- 2 ® 

^ ~ X ^ 4 V bound vortex sheet distributed over the airfoil 

chord of strength T(x-^,t). 

The circulatory velocity potential on the upper semicircle due to 

the vortex pair P , -P is 

o’ o 


<i) ie,t) 

c 


I]o ^^„~^^ /^-b)(l+cos ^ 
” f (e*+b)(l-cos o) 


and the corresponding pressure difference, from is 

P„UT^[|*+b cos 0] 

<t b sin 0 


(B.6) 


(B.7- 


p^.(o.t) = 



Tlip InnK^'nllul vrlcH'lly n I. tho circle' iiuJuccfl by I be vnrlf'x p.i Ir 


i. s 


/) 


r 


:n:b 



k , bs ii 

\ ~{p 

«) 

1 

b 

1 (r,) - ,vb cos '1 


(n.K: 


Koj’ nrblli’nry inolions^ thcT’c wil l bo n clistribulod wriko vortex 
shoot of stronKth 7^(5, t) and tho olMcct of the shed wake is obtaiiu'd 
by roplncinr: 1’^^ by 5 ; t) and inlegratinp: over the wake. l*'or 

airl'oil motion at uniform velocity starting at t = 0 tho tangential 
velocity at tho trail ing-odge is 

>HUt^ ^ 


nnd the pressure difference on the nirfoil is 


b+Utr 


Pc' 




PU 

jtbsinO 




I. ^ ^ / 

l“-b (1-c 


/f' ■x-l-b'’ 


cosO 


7 (^l^t)dP^ 


(B.IO) 


The velocity nt the IrnilinK-edge is given by (B.3) and (B.9) and the 
Kutta condition is enforced by requiring that this velocity be zero. 


2 

f _JL_ .. 


b+Ut 

1 

J 

\ 

J COS 1> - 1 

o 

^ jib 

) V 

b 


\(V^t)dV’<• 


(B.U) 


Kcluation (B.ll) relates the known downwash^ w to the unknown wake 
vortex strength 7^('’<-,t). The first integral in (B.ll) may bi’ evalualed 
^ is spcciflr*d. Thoodorson df( itu'd ono-h.ilf Uiis inl(‘gra] as 


1 .•)- 
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Q 


- - f 

7T J 


7 T w^- sin 
a 

o cos (\i - 1 


(B.12) 


giving 


Q = Ua + iib + b(i - a)5: 


^ i ’’loP ^ h V * ; If 


for the section of Fig, A- 1 , and 


Q = Ua + lib + b(J-a)a 


+ ~T ^--~T P + -(T (d)-7r)T 

TT 10^ 2tt Tt'' 10^ '' 


(B.IB) 


+ ^ (T^l(d) - ir(l - 2c)7 


(B.l 4 ) 


for the section of Fig, II-l, The T-f unctions were evaluated by 
Theodorsen [ll] and Theodorsen and Garrick [51] tabulated in 

App, C, Thus the Kutta condition^ (B,ll)^ may be written as 


Q . 


1 

2 JTb 


X* 


(B.15) 


Integrating the pressure difference (B.IO) over the chords the circu- 
latory lift and pitching moment are 

r o ^ 

’w 


P = PU 


^ - b^‘ 7^(^*,t)de* 


(B.16) 


m“ = PUb 


i* 1 


\(t^t)d?,X- . (b. 17 ) 


- 160 - 


K^ns (1^,15) and (1^.17) show tliiiL l\w circulatorv loads aro rolala d 
Lo the wako vorlt^x stron^th by Uio two iiito^rnl express Ions 


J V^. 

and 

I 2 = ^ y^^(5-st)dpx- . (B.19) 

b 

f" r 

(The hinge moments M and M' may be expressed in terms of and 

also.) Theodorsen noted that if it were assumed that the airfoil had 
undergone simple harmonic oscillations for an indefinitely long period 
then 


w*(x*,t) = w*(x> (B.20) 

U cl 

and 

7 (!•»:-. t) = 7 ^ y (B21) 

since the wake is assumed to drift downstream at the freestrcam velocity 
as shown in IT-y. The reduced frequency, k - ub/U introduced in 

(B.21) serves to indicate the relative ’unsteadiness' of the flow. T!ie 
unknown wake vortex strength may now be factored out of the in- 

tegrals and (B.15) to (B.17) become 


7 

w 
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(B.22) 



THE AIRFOIL 



( 15 . 23 ) 


(B.24) 


where F = E*/b and a = Ut/b are nondi mens ionali zed distance and time, 
respectively. 

The integrals may be identified as modified Bessel functions of 
the third kind, K^,(s), [Ref. 57, p. 22] from the integral definition 
[Ref. 58, Eq. 9.6.23] 




r(iU|) 

r(v+i) 


03 

1 


2 . V-? 


(t -I)" =^dt; 


Refv) > -i, 
Rc(s) > 0, 


(B.25) 


where T(«) if the Gamma function [57], For v = 0 


K„( = ) 


^ -si 
o 


(B.26) 


and sinci' Kj^(s) -K]^(s) 
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Therefore 


K^(s). + K^(s) = ^ 03.28) 


and the ratio of integrals in (B.23) and (B.24), defined as the Theodor- 
sen function C(ik), is 


K (ik) 

C(ik) = i 

K^(ik) + K^(ik) 




(B.29) 


( 2 ) 

The Hankel functions are given by W 

(-n/2) (k) [Ref. 58, Eq. 9,6,4]. Theodorsen [11] did not 

mention the violat ion of the condition Re (s) strictly greater than 
zero, in the application of (B.26) to (B.23) and (B,24). 

The loads acting on the airfoil may be calculated from (B.5), 

(B.14), (B.23), (B.24), and (B.29) and from similar equations for the 

hinge moments. The integrals required were evaluated by Theodorsen [11] 
and Theodorsen and Garrick [51] and are tabulated in App. C. lor the 
wing-aileron-tab section of Fig. A-1 , the loads may be writte n as 

L = L + Pb'^M X + Pb^UB x + Pb^U^K X (B.30) 

c nc" nc no- 

where the matrices giving the Mionc i rculatory ’ loads are 
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— Y 
T 10 
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V 18 


- T,„(d) 
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0 0 
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0 0 - -T 


1 H 
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i '9 




[ind the ’circulatory’ loads arc given by 


L = Pb^UC(ik)RQ 
c 


where 
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The factor Q, (n.l.'t), may be urilLoii 


Q = US X 1- bS X 


(13. :52) 


where 


\ [«' i'V 


[l, (5-»), 2, ’’ll’ 2ir'‘u'‘'^ 


The loads for the section with leading:- and i ra i 1 i nn-efine control 
surfaces (Fi(^, II-l).-may be obtained from (B.30) and (11. SI) using 
(A. 17) and A. 18). They are 


L = L ! Pb^M X t Pb^UB x + Pb^U^K x 

... - n r*“ n f ’ ~ 


(B.SS) 


who re 


M ^ V M V 

nc nc 

B - v'^B V 

iu‘ nc 

K V^K V 

nc nc 


/'ll 'tK ' Ik IvSii S Vn • S V.x] 


(B.S'l) 


II. should l)c nolid t liat not all iatiLors conipi - 1 s i n/; Mic circulalorv 
lads arc multiplied l.\ i (lU). I'his results in i-erta i n cancel lat i oiis , fVi 

-I 


I 


1 


of torms botwcen tho circulatory and noncirculatory lojuU and the sub- 
scriptH in (l^.dO) find (B.ISM) aro only doscriptivo . 



Appondix C 


aNSTh;AJ)V AKKOUYNAMIC LOADS IN TWO^DIMENSIONAl, 
INCQMPllKSSIBLK FLOW 


The aerodynamic loads acting on the section of Fifr, A-l were cal- 
culated by Theodorsen [U], and Theodorsen and Garrick rr>l] and Involve 
the follow i nr: expressions . 
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Appendix I) 


AI.THUNATIV I'; DKUIVATIONS Tlffi r.KNKUAI.r/.KI) 
TlUjODOHSKN FUNCTION 


1. W.p. Jones [20], using the concepts of bound and I'l-ee vortlclty, 
was able to show that the functions involved in Theodorsen's problem 
satisfied the modified Bessel equation. He thus avoided the restriction, 
Re(s) > 0, involved in the integral j-epresentations of tlie Bessel func- 
tions and derived the generalized Theodorsen function 


K^(s) + K^(s) 

In attempting to compare this form of C(s) with that given-by Theo- 
dorsen, Jones used the relation 


K^(s) = J (D.2) 

without regard to the restriction -TT < arg s s ^ TT. As a result, 

C(s) was evaluated using different branches of the functions Il^^\s) and 

(2) - o 

H^^ (s) in the first and second quadrants of the s-plane and Jones con- 
cluded, incorrectly, that C(s) was discontinuous across tlie i.. axis, 

2. The convolution integral may bo used to verify that the gen- 
eralized Theodorsen function is indeed the correct operator relating 
the downwash, w(s), to the induced airloads for stable airfoil motions. 
The lift due to circulation is given by 

P(«) = 2'iPbUC(s)w(s) (D.B) 

lor some o|ieiato)’, (’(s). For the' assuim.'fl form of ( 2 .BB), 
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Kj(s) 


1 K^(h) 

nnd tlio pai^ticuliu’ (lam|Hul airfoil motion 




\v( co£> wt* « Ilej^w^o^ cT-i i ) t ^ ^ 


0 ()).d) 


the? inverse Laplace transffjrm of 3 ^(h) is 


p(t') 

2TTPbU 
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27Ti 


o^+i«> 
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C(s)w(s)e^^ ds . 




The transform w(s) is }>iveu by the real part of 


(».5) 


\v(s) 


js+o) + 1^ 1 

i 2 ' 

,(si a) Tto J 


(D.6) 


Since 0 > 0, may be set equal to x.ero and with the substitution 

s ;= ik 


P(t' ) 

2TTPbU 





+ ikt» , 
e dk. 




(IJ.7) 


i'he s ymbol under the i ntep;ral Implies that the path of intoxication 
must pass below the branch pt^int at tlie origin, Ciarrick [70] showed 
that tile lift could be calculated inv arbitrary motions usinx the convo- 
lution Intexral 


✓ 



- 1 72- 


I 


t' 


" »(0)k,(t.) H. j kj(t') flililii <I,' 


(n.H) 


(If 


Alsu, ^Uk)/k nro rolat(Ml ! 70 j by 


k ( 1 1 ) - _i_ f £ 
’ 2 «i J 


(k) Ikf 


ell 


( 0 . 9 ) 


-ir> 


where the fact that k^(t’) 0 for t^ < 0 has been used* Substituting 

the expressions Kiven by (D,6) and (0,9) into (D.8) yields 
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The expression in brackets inside the last integral represents k^(t^) 
which is zero for t* < 0 and the integrals may be interchanged, giving 
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Then (D.lo) is 
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UlSCUSSION OF TIIF GENFRAI.IZKI) THEODOUSKN FUNCTION AND 
UNSTEADY AERODYNAMICS FOR ARBITRARY MOTIONS 


The study of unsteady airloads due to transient motions was pioneered 
by Wagner [ S] who calculated the lift on an airfoil started impulsively 
from rest. The resulting lift function kj^(t') is known as Wagner's func- 
tion and has not been successfully evaluated in terms of elementary func- 
tions. Due to the linearity of the governing partial differential equa- 
tions, it was recognized that superposition of elementary solutions could 
be used to calculate unsteady airloads for arbitrary motions. Garrick 
[70] used the convolution Integral to write the lift due to motion w(t' ) 
as 


P(f ) 
2nPbU 


t« 

w(0)k^(t') + J k^(t'-t|) 


dw(-t|) 




dt| 


(E.l) 


Garrick [701 also showed that k^ (t?) and C(ik)/k were a Fourier trans- 
form pair, 


00 

C(ik) = ik f k^(t')c dt' 


(E.2) 
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In (K,2) Ihi? fad that k^(LM ^ H i’or t* < 0 has bcuMi uschI and the path 
nl i nl(*nra t ion in (E,.M) must pass bi'low I lie:' sinji;ul ar 11 y ol‘ t, he inte^i’and 
at k 0. 
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Tiu‘ .appl icaLion Laplaco trans form tocjlmiqucs to uiislfsuly atoH)- 
(lynaiiilc intc^l^ral ociuations wn.s si 4 'p;ost:ocl by U/l\ Jonos [29] and Seaj’s [bb] 
usod the LoclmiciiU' to obtain new solutions to Wagner’s problem (indlcial 
lift diu' to impulsive plunp;ing), Kussner's |)roblcm (indiclal lift, 
due to pc'iu'Lration of a sharp~ed 4 ;ed ^ust), and thc^ oscillating airfoil 
})roblem, wSears * presentation is essentially a derivation of the general- 
i/ed Theodorsen function although this aspect is not discussed in Ref, bb 
and was apparently not pursued. 

It is interesting to note that the early references in the field do 
not mention the restrictions on the existence of the integrals upon which 
the theory is based, Sohngen [88] was apparently the first to recognize 
the effect of the branch cut of C(s) upon the loads. He noted that 

diverging airfoil motions led to airloads which behaved asymptotically 

st * 

as e' while converging airfoil motions led to asymptotic loads pro- 
portional to l/t\ These correspond to the rational and nonrational 
])ortions identified in the text. This difference was bothersome and it 
appeared to correlate with the restriction upon the existence of the 
integrals in question (viz,, Re(s) > 0), leading to the conclusion that 
C(s) could not be extended into the left half-plane. 

The first attempt to evaluate the Theodorsen function for complex 
values of k was by W.P, Jones [20] who concluded, incorrectly, that 
C(s) was discontinuous across the imaginary axis. Thus, he concluded 
that C(ik) could be generalized for divergent oscillations (Fig, II-3a) 
but was invalid for convergent oscillations (Fig. II-3b), This reasoning 
was reinforced by the fact that Theodorsen had been forced to ‘assume an 
explicit form for the airfoil motion and wake vortex distribution (B,20, 
13.21) in order to obtain a solution. This fact may be the source of the 
Confusion wherein the Theodorsen function is interpreted as a time domain 
operator rather than a frequency domain operator (e.g,, Rc’f. 22), 

During this period, calculations were made of unsteady loads using 
the convolution integral (K,l) with the indicial function approximated 
by sums of exponential time factoi’s as shc^wn by R,T, Jones [29], [B]. 

The exponential appi'ox imations were cai)able of being evaluated for 




— 1 7 


I 


ai’l)itrury motion, luul Goliunl and IaiKo [30) published root looi ,»J aero- 

elastic modes, Tlie Laplace transform of the exponontinJ approximation 

to interprotcd ns an ad hoc ^oncrali^od Thoodorson function 

and comparisons of numerical calculations [25] usin^ such functions and 

the exact tabu^lated k (t*) function led Luke and DciiKlcr [21] to the 
i f 

conclusion that G(ik) could be extended to the entire s-plane. However, 
their arf^umont based upon analytic continuation, did not seem convincing 
in light of the abdVe discussion and it was rejected in a series of 
articles [22], [23]^' [24], [25], and [26], At the heart of the discussion 
was the requirement in Theodorsen’s derivation of assuming an explicit 
form of airfoil motion and wake vortex distribution (i,e,, oscillatory 
divergent and infinite extent) in order to evaluate the resulting. inte~ 
grals. It seemed contrary to reason to claim that the resulting function 
C(ik) was valid for damped motions when the derivation of the function 
required just the opposite assumption. Of course, the assumption of an 
explicit motion is not required and the derivation of the generalized 
Theodorsen function using Laplace transform techniques is given in the 
text, 

It would appear that the difficulties with the generalized Theodorsen 
function influenced the subsequent development of compressible finite-wing 
aerodynamic load calculations. These techniques [15], [17], [62] and [89] 
invariably begin with the assumption of simple harmonic oscillations, 
although the text shows that this assumption is not necessary. It must 
be recognized that oscillatory loads are entirely adequate forthe estab- 
lishment of flutter boundaries and, until the advent of active aeroclastic 
control schemes, there was little requirement for loads due to arbitrary 
motions , 

Morino [36] has derived a new formulation of the unsteady aerodynamic 
loading problem based upon the Green function solution and claims that 
this formulation is the only technique capable of analyzing loads due 
to arbitrary motions. Presumably -this claim is based upon the fac4;— that 
Morino’ s theory analyzes finite thickness wings and does not encounter 
the singularities inherent in flat plate theories. However, the text 
of this thesis shows that these singularities do not restrict consideration 
to oscillatory motions and Morino’s claim is false, 
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Appendix F 

■SERIES EXPANSIONS OF BESSEL FUNCTIONS 


The following series expansions are given In Ch. 9 of Abramowitz and 
Ste^^un [58], 

1. Ascending power series 
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